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

IndexAndRefine: Another iteration

Browse Source
This commit is contained in:
leonarski_f
2025-11-30 13:20:26 +01:00
parent 69afeb488f
commit cb983f5c74
5 changed files with 280 additions and 196 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/CMakeLists.txt
+3 -1
View File
@@ -6,7 +6,9 @@ ADD_LIBRARY(JFJochImageAnalysis STATIC
IndexAndRefine.cpp
IndexAndRefine.h
dark_mask_analysis/DarkMaskAnalysis.cpp
dark_mask_analysis/DarkMaskAnalysis.h)
dark_mask_analysis/DarkMaskAnalysis.h
RotationIndexer.cpp
RotationIndexer.h)
FIND_PACKAGE(Eigen3 3.4 REQUIRED NO_MODULE) # provides Eigen3::Eigen
image_analysis/IndexAndRefine.cpp
+45 -163
View File
@@ -11,99 +11,17 @@
#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()),
: 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;
}
}
}
if (indexer && x.GetGoniometer() && x.GetIndexingSettings().GetRotationIndexing())
rotation_indexer = std::make_unique<RotationIndexer>(x, *indexer);
}
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;
}
}
if (rotation_indexer)
rotation_indexer->SetLattice(lattice);
}
void IndexAndRefine::ProcessImage(DataMessage &msg,
@@ -125,13 +43,24 @@ void IndexAndRefine::ProcessImage(DataMessage &msg,
bool beam_center_updated = false;
if (!axis_) {
if (rotation_indexer) {
auto result = rotation_indexer->ProcessImage(msg.number, msg.spots);
if (result) {
lattice_candidate = result->lattice;
experiment_copy.BeamX_pxl(result->geom.GetBeamX_pxl())
.BeamY_pxl(result->geom.GetBeamY_pxl())
.DetectorDistance_mm(result->geom.GetDetectorDistance_mm());
symmetry.centering = result->search_result.centering;
symmetry.niggli_class = result->search_result.niggli_class;
symmetry.crystal_system = result->search_result.system;
}
} else {
// 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_));
recip.push_back(i.ReciprocalCoord(geom_));
}
auto indexer_result = indexer_->Run(experiment, recip).get();
@@ -163,72 +92,40 @@ void IndexAndRefine::ProcessImage(DataMessage &msg,
latt = sym_result.conventional;
}
XtalOptimizerData data{
.geom = geom_,
.latt = latt,
.crystal_system = symmetry.crystal_system,
.min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots(),
.refine_beam_center = true,
.refine_distance_mm = false
};
lattice_candidate = latt;
}
} else {
std::unique_lock ul(m);
const auto rot = axis_->GetTransformation(msg.number);
if (symmetry.crystal_system == gemmi::CrystalSystem::Trigonal)
data.crystal_system = gemmi::CrystalSystem::Hexagonal;
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_));
}
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;
}
break;
}
accumulated_images++;
last_accumulated_image = msg.number;
lattice_candidate = data.latt;
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;
experiment_copy.BeamX_pxl(updated_geom_.GetBeamX_pxl())
.BeamY_pxl(updated_geom_.GetBeamY_pxl())
.DetectorDistance_mm(updated_geom_.GetDetectorDistance_mm());
if (beam_center_updated) {
msg.beam_corr_x = data.beam_corr_x;
msg.beam_corr_y = data.beam_corr_y;
}
}
}
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 (AnalyzeIndexing(msg, experiment, data.latt)) {
if (AnalyzeIndexing(msg, experiment, *lattice_candidate)) {
msg.lattice_type = symmetry;
float ewald_dist_cutoff = 0.001f;
@@ -239,11 +136,6 @@ void IndexAndRefine::ProcessImage(DataMessage &msg,
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);
@@ -262,17 +154,7 @@ void IndexAndRefine::ProcessImage(DataMessage &msg,
}
std::optional<RotationIndexerResult> IndexAndRefine::Finalize(bool retry) {
if (axis_) {
std::unique_lock ul(m);
if (!indexed_lattice || retry)
TryIndexRotationData();
if (indexed_lattice)
return RotationIndexerResult{
.lattice = indexed_lattice.value(),
.search_result = search_result_,
.geom = geom_
};
}
if (rotation_indexer)
return rotation_indexer->Finalize(retry);
return {};
}
image_analysis/IndexAndRefine.h
+5 -32
View File
@@ -1,8 +1,8 @@
// 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
#ifndef JFJOCH_INDEXANDREFINE_H
#define JFJOCH_INDEXANDREFINE_H
#include <vector>
#include <mutex>
@@ -12,46 +12,19 @@
#include "bragg_integration/BraggPrediction.h"
#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;
DiffractionGeometry geom;
};
#include "RotationIndexer.h"
class IndexAndRefine {
constexpr static int64_t max_spots = 60000;
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_;
DiffractionGeometry updated_geom_;
LatticeSearchResult search_result_;
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 = 1;
int64_t first_image_to_try_indexing = INT64_MAX;
void TryIndexRotationData();
std::unique_ptr<RotationIndexer> rotation_indexer;
public:
IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool *indexer);
void SetLattice(const CrystalLattice &lattice);
@@ -60,4 +33,4 @@ public:
};
#endif //JFJOCH_SPOTROTATIONSTORE_H
#endif //JFJOCH_INDEXANDREFINE_H
image_analysis/RotationIndexer.cpp
+163
View File
@@ -0,0 +1,163 @@
// 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/RotationIndexer.h
+64
View File
@@ -0,0 +1,64 @@
// 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
#include <vector>
#include <mutex>
#include "../common/DiffractionSpot.h"
#include "../common/DiffractionExperiment.h"
#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;
DiffractionGeometry geom;
};
class RotationIndexer {
mutable std::mutex m;
const DiffractionExperiment& experiment;
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 bool index_ice_rings;
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;
int64_t image_stride;
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<RotationIndexerResult> ProcessImage(int64_t image, const std::vector<SpotToSave>& spots);
std::optional<RotationIndexerResult> Finalize(bool retry);
};
#endif //JFJOCH_ROTATIONINDEXER_H