diff --git a/image_analysis/RotationIndexer.cpp b/image_analysis/RotationIndexer.cpp index ba0ad0a6..f3ff9372 100644 --- a/image_analysis/RotationIndexer.cpp +++ b/image_analysis/RotationIndexer.cpp @@ -6,6 +6,7 @@ #include "geom_refinement/XtalOptimizer.h" #include "indexing/FFTIndexer.h" #include "lattice_search/LatticeSearch.h" +#include RotationIndexer::RotationIndexer(const DiffractionExperiment &x, IndexerThreadPool &indexer) : experiment(x), @@ -15,7 +16,7 @@ RotationIndexer::RotationIndexer(const DiffractionExperiment &x, IndexerThreadPo updated_geom_(geom_), indexer_(indexer) { if (axis_) { - float angle_norm_deg = std::fabs(axis_->GetIncrement_deg()); + angle_norm_deg = std::fabs(axis_->GetIncrement_deg()); if (angle_norm_deg < 1e-6) { // Guard against rotation close to zero axis_ = std::nullopt; @@ -97,6 +98,12 @@ void RotationIndexer::TryIndex() { indexed_lattice = data.latt; updated_geom_ = data.geom; } + PredictionRotationSettings pred_settings{ + .high_res_A = 1.0, + .max_hkl = 100, + .centering = search_result_.centering + }; + predictions = PredictRotationHKLs(experiment, data.latt, pred_settings); } } @@ -144,3 +151,7 @@ std::optional RotationIndexer::GetLattice() { .geom = updated_geom_ }; } + +const std::vector &RotationIndexer::GetPredictions() const { + return predictions; +} diff --git a/image_analysis/RotationIndexer.h b/image_analysis/RotationIndexer.h index be6fde20..14b2efd5 100644 --- a/image_analysis/RotationIndexer.h +++ b/image_analysis/RotationIndexer.h @@ -11,6 +11,7 @@ #include "../common/DiffractionExperiment.h" #include "indexing/IndexerThreadPool.h" #include "lattice_search/LatticeSearch.h" +#include "bragg_prediction/BraggPredictionRotation.h" // RotationIndexer works as following: // 1. First accumulates spot results from rotation images (only images within a certain stride are included) @@ -33,6 +34,7 @@ class RotationIndexer { const bool index_ice_rings; bool indexing_tried = false; + float angle_norm_deg = 1.0f; std::vector v_; std::vector coords_; @@ -51,11 +53,13 @@ class RotationIndexer { std::optional indexed_lattice; + std::vector predictions; void TryIndex(); public: RotationIndexer(const DiffractionExperiment& x, IndexerThreadPool& indexer); std::optional ProcessImage(int64_t image, const std::vector& spots); std::optional GetLattice(); + const std::vector &GetPredictions() const; }; #endif //JFJOCH_ROTATIONINDEXER_H \ No newline at end of file diff --git a/image_analysis/bragg_prediction/BraggPredictionRotation.cpp b/image_analysis/bragg_prediction/BraggPredictionRotation.cpp index 5c55ded4..8ad41404 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRotation.cpp +++ b/image_analysis/bragg_prediction/BraggPredictionRotation.cpp @@ -11,66 +11,11 @@ #include "../../common/JFJochException.h" #include "../bragg_integration/SystematicAbsence.h" -BraggPredictionRotation::BraggPredictionRotation(int max_reflections) - : BraggPrediction(max_reflections) {} +std::vector PredictRotationHKLs(const DiffractionExperiment &experiment, + const CrystalLattice &lattice, + const PredictionRotationSettings &settings) { + std::vector ret; -namespace { - inline float deg_to_rad(float deg) { - return deg * (static_cast(M_PI) / 180.0f); - } - inline float rad_to_deg(float rad) { - return rad * (180.0f / static_cast(M_PI)); - } - - // Solve A cos(phi) + B sin(phi) + D = 0, return solutions in [phi0, phi1] - // Returns 0..2 solutions. - inline int solve_trig(float A, float B, float D, - float phi0, float phi1, - float out_phi[2]) { - if (phi1 < phi0) std::swap(phi0, phi1); - - const float R = std::sqrt(A*A + B*B); - if (!(R > 0.0f)) - return 0; - - const float rhs = -D / R; - if (rhs < -1.0f || rhs > 1.0f) - return 0; - - const float phi_ref = std::atan2(B, A); - const float delta = std::acos(rhs); - - float s1 = phi_ref + delta; - float s2 = phi_ref - delta; - - // Shift candidates by +/- 2*pi so they land near the interval. - const float two_pi = 2.0f * static_cast(M_PI); - auto shift_near = [&](float x) { - while (x < phi0 - two_pi) x += two_pi; - while (x > phi1 + two_pi) x -= two_pi; - return x; - }; - - s1 = shift_near(s1); - s2 = shift_near(s2); - - int n = 0; - if (s1 >= phi0 && s1 <= phi1) out_phi[n++] = s1; - if (s2 >= phi0 && s2 <= phi1) { - if (n == 0 || std::fabs(s2 - out_phi[0]) > 1e-6f) out_phi[n++] = s2; - } - return n; - } - - // Convert angle to fractional image_number using goniometer start/increment - inline float phi_deg_to_image_number(float phi_deg, float start_deg, float increment_deg) { - return (phi_deg - start_deg) / increment_deg; - } -} // namespace - -int BraggPredictionRotation::Calc(const DiffractionExperiment &experiment, - const CrystalLattice &lattice, - const BraggPredictionSettings &settings) { const auto gon_opt = experiment.GetGoniometer(); if (!gon_opt.has_value()) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, @@ -78,15 +23,12 @@ int BraggPredictionRotation::Calc(const DiffractionExperiment &experiment, const GoniometerAxis& gon = *gon_opt; const auto geom = experiment.GetDiffractionGeometry(); - const float det_width_pxl = static_cast(experiment.GetXPixelsNum()); - const float det_height_pxl = static_cast(experiment.GetYPixelsNum()); const float one_over_dmax = 1.0f / settings.high_res_A; const float one_over_dmax_sq = one_over_dmax * one_over_dmax; const float wavelength = geom.GetWavelength_A(); const float one_over_wavelength = 1.0f / wavelength; - const float one_over_wavelength_sq = one_over_wavelength * one_over_wavelength; const Coord Astar = lattice.Astar(); const Coord Bstar = lattice.Bstar(); @@ -100,128 +42,96 @@ int BraggPredictionRotation::Calc(const DiffractionExperiment &experiment, const float beam_y = geom.GetBeamY_pxl(); const float det_distance = geom.GetDetectorDistance_mm(); const float pixel_size = geom.GetPixelSize_mm(); - const float coeff_const = det_distance / pixel_size; - - // Determine prediction interval in spindle angle. - const float wedge_deg = gon.GetWedge_deg(); - - const float phi0_deg = gon.GetAngle_deg(static_cast(settings.image_number)); - const float phi1_deg = phi0_deg + wedge_deg; - - const float phi0 = deg_to_rad(phi0_deg); - const float phi1 = deg_to_rad(phi1_deg); - - // Convert mosaicity to radians (assumes settings has this field, or defaults to small value) - const float mosaicity_rad = deg_to_rad(settings.mosaicity_deg); - const float half_mosaicity = mosaicity_rad * 0.5f; + const float F = det_distance / pixel_size; const Coord m2 = gon.GetAxis().Normalize(); + const Coord m1 = (m2 % S0).Normalize(); + const Coord m3 = (m1 % m2).Normalize(); + + const float m2_S0 = m2 * S0; + const float m3_S0 = m3 * S0; int i = 0; - for (int h = -settings.max_hkl; h <= settings.max_hkl; ++h) { + for (int32_t h = -settings.max_hkl; h <= settings.max_hkl; ++h) { const float Ah_x = Astar.x * h; const float Ah_y = Astar.y * h; const float Ah_z = Astar.z * h; - for (int k = -settings.max_hkl; k <= settings.max_hkl; ++k) { + for (int32_t k = -settings.max_hkl; k <= settings.max_hkl; ++k) { const float AhBk_x = Ah_x + Bstar.x * k; const float AhBk_y = Ah_y + Bstar.y * k; const float AhBk_z = Ah_z + Bstar.z * k; - for (int l = -settings.max_hkl; l <= settings.max_hkl; ++l) { + for (int32_t l = -settings.max_hkl; l <= settings.max_hkl; ++l) { if (systematic_absence(h, k, l, settings.centering)) continue; - if (i >= max_reflections) - continue; - const float p0_x = AhBk_x + Cstar.x * l; const float p0_y = AhBk_y + Cstar.y * l; const float p0_z = AhBk_z + Cstar.z * l; - const float p0_sq = p0_x * p0_x + p0_y * p0_y + p0_z * p0_z; + const Coord p0{p0_x, p0_y, p0_z}; + + const float p0_sq = p0 * p0; if (p0_sq <= 0.0f || p0_sq > one_over_dmax_sq) continue; - const float p0_par_s = p0_x * m2.x + p0_y * m2.y + p0_z * m2.z; - const float p0_par_x = m2.x * p0_par_s; - const float p0_par_y = m2.y * p0_par_s; - const float p0_par_z = m2.z * p0_par_s; + const float p0_m1 = p0 * m1; + const float p0_m2 = p0 * m2; + const float p0_m3 = p0 * m3; - const float p0_perp_x = p0_x - p0_par_x; - const float p0_perp_y = p0_y - p0_par_y; - const float p0_perp_z = p0_z - p0_par_z; + const float rho_sq = p0_sq - (p0_m2 * p0_m2); - const float p0_perp_sq = p0_perp_x * p0_perp_x + p0_perp_y * p0_perp_y + p0_perp_z * p0_perp_z; - if (p0_perp_sq < 1e-12f) + const float p_m3 = (- p0_sq / 2 - p0_m2 * m2_S0) / m3_S0; + const float p_m2 = p0_m2; + const float p_m1_opt[2] = { + std::sqrt(rho_sq - p_m3 * p_m3), + -std::sqrt(rho_sq - p_m3 * p_m3) + }; + + // No solution for Laue equations + if ((rho_sq < p_m3 * p_m3) || (p0_sq > 4 * S0 * S0)) continue; - const float p_x = S0.x + p0_par_x; - const float p_y = S0.y + p0_par_y; - const float p_z = S0.z + p0_par_z; + for (const auto& p_m1 : p_m1_opt) { + const float cosphi = (p_m1 * p0_m1 + p_m3 * p0_m3) / rho_sq; + const float sinphi = (p_m1 * p0_m3 - p_m3 * p0_m1) / rho_sq; + Coord S = S0 + m1 * p_m1 + m2 * p_m2 + m3 * p_m3; - const float m2_x_gperp_x = m2.y * p0_perp_z - m2.z * p0_perp_y; - const float m2_x_gperp_y = m2.z * p0_perp_x - m2.x * p0_perp_z; - const float m2_x_gperp_z = m2.x * p0_perp_y - m2.y * p0_perp_x; + float phi = -1.0f * std::atan2(sinphi, cosphi) * 180.0f / M_PI; + if (phi < 0.0f) phi += 360.0f; + const float lorentz_reciprocal = std::fabs(m2 * (S % S0))/(S*S0); - // Trig equation coefficients for solving the Laue condition in φ: - // A cosφ + B sinφ + D = 0 - const float A_trig = 2.0f * (p_x * p0_perp_x + p_y * p0_perp_y + p_z * p0_perp_z); - const float B_trig = 2.0f * (p_x * m2_x_gperp_x + p_y * m2_x_gperp_y + p_z * m2_x_gperp_z); - const float D_trig = (p_x * p_x + p_y * p_y + p_z * p_z) + p0_perp_sq - one_over_wavelength_sq; - - float sols[2]{}; - const int nsol = solve_trig(A_trig, B_trig, D_trig, phi0 - half_mosaicity, phi1 + half_mosaicity, sols); - - for (int si = 0; si < nsol; ++si) { - if (i >= max_reflections) break; - - const float phi = sols[si]; - - const float cos_p = cosf(phi); - const float sin_p = sinf(phi); - - // Rodrigues' rotation formula: g = g_par + cos(phi)*g_perp + sin(phi)*(w x g_perp) - const float p_star_x = p0_par_x + cos_p * p0_perp_x + sin_p * m2_x_gperp_x; - const float p_star_y = p0_par_y + cos_p * p0_perp_y + sin_p * m2_x_gperp_y; - const float p_star_z = p0_par_z + cos_p * p0_perp_z + sin_p * m2_x_gperp_z; - - const float S_x = S0.x + p_star_x; - const float S_y = S0.y + p_star_y; - const float S_z = S0.z + p_star_z; - - // Inlined RecipToDector - const float S_rot_x = rot[0] * S_x + rot[1] * S_y + rot[2] * S_z; - const float S_rot_y = rot[3] * S_x + rot[4] * S_y + rot[5] * S_z; - const float S_rot_z = rot[6] * S_x + rot[7] * S_y + rot[8] * S_z; + // Inlined RecipToDector with rot1 and rot2 (rot3 = 0) + // Apply rotation matrix transpose + float S_rot_x = rot[0] * S.x + rot[1] * S.y + rot[2] * S.z; + float S_rot_y = rot[3] * S.x + rot[4] * S.y + rot[5] * S.z; + float S_rot_z = rot[6] * S.x + rot[7] * S.y + rot[8] * S.z; if (S_rot_z <= 0) continue; - const float coeff = coeff_const / S_rot_z; - const float x = beam_x + S_rot_x * coeff; - const float y = beam_y + S_rot_y * coeff; + float x = beam_x + F * S_rot_x / S_rot_z; + float y = beam_y + F * S_rot_y / S_rot_z; - if (x < 0.0f || x >= det_width_pxl || y < 0.0f || y >= det_height_pxl) - continue; - float dist_ewald_sphere = std::fabs(sqrtf(S_x * S_x + S_y * S_y + S_z * S_z) - one_over_wavelength); - float d = 1.0f / sqrtf(p0_sq); - const float phi_deg = rad_to_deg(phi); - reflections[i] = Reflection{ + ret.emplace_back(PredictionRotationResult{ + .angle_deg = phi, + .lorentz_reciprocal = lorentz_reciprocal, .h = h, .k = k, .l = l, - .angle_deg = phi_deg, - .predicted_x = x, - .predicted_y = y, - .d = d, - .dist_ewald = dist_ewald_sphere - }; - ++i; + .x = x, + .y = y + }); } } } } - return i; -} \ No newline at end of file + + std::ranges::sort(ret, + [](const PredictionRotationResult& a, const PredictionRotationResult& b) { + return a.angle_deg < b.angle_deg; + }); + return ret; +} diff --git a/image_analysis/bragg_prediction/BraggPredictionRotation.h b/image_analysis/bragg_prediction/BraggPredictionRotation.h index a411f3d5..40b9a607 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRotation.h +++ b/image_analysis/bragg_prediction/BraggPredictionRotation.h @@ -12,15 +12,22 @@ #include "../../common/Reflection.h" #include "BraggPrediction.h" - -class BraggPredictionRotation : public BraggPrediction { -public: - BraggPredictionRotation(int max_reflections = 200*200*200); - - int Calc(const DiffractionExperiment &experiment, - const CrystalLattice &lattice, - const BraggPredictionSettings &settings) override; +struct PredictionRotationSettings { + float high_res_A; + int32_t max_hkl = 100; + char centering = 'P'; }; +struct PredictionRotationResult { + float angle_deg; + float lorentz_reciprocal; + int32_t h,k,l; + float x,y; +}; + +std::vector PredictRotationHKLs(const DiffractionExperiment &experiment, + const CrystalLattice &lattice, + const PredictionRotationSettings &settings); + #endif //JFJOCH_BRAGGPREDICTIONROTATION_H \ No newline at end of file