diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index dc8094d3..8bc8b50a 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -188,7 +188,8 @@ void IndexAndRefine::ProcessImage(DataMessage &msg, if (!outcome.lattice_candidate) return; - if (!AnalyzeIndexing(msg, outcome.experiment, *outcome.lattice_candidate, outcome.experiment.GetGoniometer())) + if (!AnalyzeIndexing(msg, outcome.experiment, *outcome.lattice_candidate, + outcome.experiment.GetGoniometer(), rotation_indexer.get() != nullptr)) return; msg.lattice_type = outcome.symmetry; diff --git a/image_analysis/indexing/AnalyzeIndexing.cpp b/image_analysis/indexing/AnalyzeIndexing.cpp index 14d95e3b..2937dcad 100644 --- a/image_analysis/indexing/AnalyzeIndexing.cpp +++ b/image_analysis/indexing/AnalyzeIndexing.cpp @@ -116,12 +116,58 @@ namespace { return rad_to_deg(best_phi); } + + std::optional CalcMosaicity(const DiffractionExperiment& experiment, + const std::vector &spots, + const Coord &astar, const Coord &bstar, const Coord &cstar) { + const auto &axis = experiment.GetGoniometer(); + if (axis.has_value()) { + const Coord w = axis->GetAxis().Normalize(); + const Coord S0 = experiment.GetScatteringVector(); + const float wedge_deg = axis->GetWedge_deg(); + const float start_deg = axis->GetStart_deg(); + const float inc_deg = axis->GetIncrement_deg(); + + double sum_sq = 0.0; + int count = 0; + + for (const auto &s: spots) { + if (!s.indexed) + continue; + + // Observed angle: use frame center + const float image_center = static_cast(s.image) + 0.5f; + const float phi_obs_deg = start_deg + inc_deg * image_center; + + // g0 at phi=0 assumption + const Coord g0 = astar * static_cast(s.h) + + bstar * static_cast(s.k) + + cstar * static_cast(s.l); + + // Local solve window: +/- 1 wedge (easy/robust first try) + const auto phi_pred_deg_opt = predict_phi_deg_local(g0, S0, w, phi_obs_deg, wedge_deg); + if (!phi_pred_deg_opt.has_value()) + continue; + + float dphi = wrap_deg_pm180(phi_obs_deg - phi_pred_deg_opt.value()); + sum_sq += static_cast(dphi) * static_cast(dphi); + count++; + } + + if (count > 0) { + return static_cast(std::sqrt(sum_sq / static_cast(count))); + } + } + return std::nullopt; + } + } // namespace bool AnalyzeIndexing(DataMessage &message, const DiffractionExperiment &experiment, const CrystalLattice &latt, - const std::optional &rotation_axis) { + const std::optional &rotation_axis, + bool rotation_indexed) { std::vector indexed_spots(message.spots.size()); // Check spots @@ -172,8 +218,9 @@ bool AnalyzeIndexing(DataMessage &message, nspots_ref++; } - if (nspots_indexed >= viable_cell_min_spots - && nspots_indexed >= std::lround(min_percentage_spots * nspots_ref)) { + if (rotation_indexed + || (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 false; @@ -187,50 +234,7 @@ bool AnalyzeIndexing(DataMessage &message, message.spot_count_indexed = nspots_indexed; message.indexing_lattice = latt; message.indexing_unit_cell = latt.GetUnitCell(); - - message.mosaicity_deg = std::nullopt; - if (rotation_axis.has_value()) { - const auto gon_opt = experiment.GetGoniometer(); - if (gon_opt.has_value()) { - const auto &gon = *gon_opt; - - const Coord w = rotation_axis->GetAxis().Normalize(); - const Coord S0 = geom.GetScatteringVector(); - const float wedge_deg = gon.GetWedge_deg(); - const float start_deg = gon.GetStart_deg(); - const float inc_deg = gon.GetIncrement_deg(); - - double sum_sq = 0.0; - int count = 0; - - for (const auto &s: message.spots) { - if (!s.indexed) - continue; - - // Observed angle: use frame center - const float image_center = static_cast(s.image) + 0.5f; - const float phi_obs_deg = start_deg + inc_deg * image_center; - - // g0 at phi=0 assumption - const Coord g0 = astar * static_cast(s.h) - + bstar * static_cast(s.k) - + cstar * static_cast(s.l); - - // Local solve window: +/- 1 wedge (easy/robust first try) - const auto phi_pred_deg_opt = predict_phi_deg_local(g0, S0, w, phi_obs_deg, wedge_deg); - if (!phi_pred_deg_opt.has_value()) - continue; - - float dphi = wrap_deg_pm180(phi_obs_deg - phi_pred_deg_opt.value()); - sum_sq += static_cast(dphi) * static_cast(dphi); - count++; - } - - if (count > 0) { - message.mosaicity_deg = static_cast(std::sqrt(sum_sq / static_cast(count))); - } - } - } + message.mosaicity_deg = CalcMosaicity(experiment, message.spots, astar, bstar, cstar); return true; } diff --git a/image_analysis/indexing/AnalyzeIndexing.h b/image_analysis/indexing/AnalyzeIndexing.h index 6757f3be..deb2b335 100644 --- a/image_analysis/indexing/AnalyzeIndexing.h +++ b/image_analysis/indexing/AnalyzeIndexing.h @@ -13,7 +13,8 @@ constexpr static float min_percentage_spots = 0.20f; bool AnalyzeIndexing(DataMessage &message, const DiffractionExperiment &experiment, const CrystalLattice &latt, - const std::optional &rotation_axis); + const std::optional &rotation_axis, + bool rotation_indexed); #endif //JFJOCH_ANALYZEINDEXING_H \ No newline at end of file