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

AnalyzeIndexing: Handle ice rings in counting percentage of indexed spots

Browse Source
This commit is contained in:
leonarski_f
2025-12-17 20:24:58 +01:00
parent 699e78e172
commit 63d5a2dc4b
1 changed files with 16 additions and 11 deletions
Download Patch File Download Diff File Expand all files Collapse all files

27
image_analysis/indexing/AnalyzeIndexing.cpp
View File

@@ -6,6 +6,7 @@
#include "AnalyzeIndexing.h"
#include "FitProfileRadius.h"
#include "spdlog/fmt/bundled/format.h"
namespace {
inline bool ok(float x) {
@@ -122,10 +123,7 @@ bool AnalyzeIndexing(DataMessage &message,
const DiffractionExperiment &experiment,
const CrystalLattice &latt,
const std::optional<GoniometerAxis> &rotation_axis) {
size_t nspots = message.spots.size();
uint64_t indexed_spot_count = 0;
std::vector<uint8_t> indexed_spots(nspots);
std::vector<uint8_t> indexed_spots(message.spots.size());
// Check spots
const Coord a = latt.Vec0();
@@ -136,10 +134,15 @@ bool AnalyzeIndexing(DataMessage &message,
const Coord bstar = latt.Bstar();
const Coord cstar = latt.Cstar();
const bool index_ice_ring = experiment.GetIndexingSettings().GetIndexIceRings();
const auto geom = experiment.GetDiffractionGeometry();
const auto indexing_tolerance = experiment.GetIndexingSettings().GetTolerance();
const auto indexing_tolerance_sq = indexing_tolerance * indexing_tolerance;
const auto viable_cell_min_spots = experiment.GetIndexingSettings().GetViableCellMinSpots();
size_t nspots_ref = 0;
size_t nspots_indexed = 0;
// identify indexed spots
for (int i = 0; i < message.spots.size(); i++) {
auto recip = message.spots[i].ReciprocalCoord(geom);
@@ -157,22 +160,24 @@ bool AnalyzeIndexing(DataMessage &message,
Coord recip_pred = std::round(h_fp) * astar + std::round(k_fp) * bstar + std::round(l_fp) * cstar;
// See indexing_peak_check() in peaks.c in CrystFEL
if (norm_sq < indexing_tolerance * indexing_tolerance) {
indexed_spot_count++;
if (norm_sq < indexing_tolerance_sq) {
if (index_ice_ring || !message.spots[i].ice_ring)
nspots_indexed++;
indexed_spots[i] = 1;
message.spots[i].dist_ewald_sphere = geom.DistFromEwaldSphere(recip_pred);
message.spots[i].h = std::lround(h_fp);
message.spots[i].k = std::lround(k_fp);
message.spots[i].l = std::lround(l_fp);
}
if (index_ice_ring || !message.spots[i].ice_ring)
nspots_ref++;
}
auto spot_count_threshold = std::max<int64_t>(viable_cell_min_spots, std::lround(min_percentage_spots * nspots));
if (indexed_spot_count >= spot_count_threshold) {
if (nspots_indexed >= viable_cell_min_spots
&& nspots_indexed >= std::lround(min_percentage_spots * nspots_ref)) {
auto uc = latt.GetUnitCell();
if (!ok(uc.a) || !ok(uc.b) || !ok(uc.c) || !ok(uc.alpha) || !ok(uc.beta) || !ok(uc.gamma))
return {};
return false;
message.indexing_result = true;
assert(indexed_spots.size() == message.spots.size());
@@ -180,7 +185,7 @@ bool AnalyzeIndexing(DataMessage &message,
message.spots[i].indexed = indexed_spots[i];
message.profile_radius = FitProfileRadius(message.spots);
message.spot_count_indexed = indexed_spot_count;
message.spot_count_indexed = nspots_indexed;
message.indexing_lattice = latt;
message.indexing_unit_cell = latt.GetUnitCell();