diff --git a/image_analysis/indexing/AnalyzeIndexing.cpp b/image_analysis/indexing/AnalyzeIndexing.cpp index e503c81f..74cdb170 100644 --- a/image_analysis/indexing/AnalyzeIndexing.cpp +++ b/image_analysis/indexing/AnalyzeIndexing.cpp @@ -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 &rotation_axis) { - size_t nspots = message.spots.size(); - - uint64_t indexed_spot_count = 0; - std::vector indexed_spots(nspots); + std::vector 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(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();