diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 1851280e..2bfd1989 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -141,9 +141,9 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, BraggPredictionSettings settings_prediction{ .high_res_A = experiment.GetBraggIntegrationSettings().GetDMinLimit_A(), - .ewald_dist_cutoff = ewald_dist_cutoff, .max_hkl = 100, - .centering = outcome.symmetry.centering + .centering = outcome.symmetry.centering, + .ewald_dist_cutoff = ewald_dist_cutoff, }; auto nrefl = prediction.Calc(experiment, *outcome.lattice_candidate, settings_prediction); diff --git a/image_analysis/bragg_prediction/BraggPrediction.cpp b/image_analysis/bragg_prediction/BraggPrediction.cpp index 3a7ae4bb..00f2389b 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.cpp +++ b/image_analysis/bragg_prediction/BraggPrediction.cpp @@ -14,19 +14,19 @@ const std::vector &BraggPrediction::GetReflections() const { int BraggPrediction::Calc(const DiffractionExperiment &experiment, const CrystalLattice &lattice, const BraggPredictionSettings &settings) { - auto geom = experiment.GetDiffractionGeometry(); - auto det_width_pxl = static_cast(experiment.GetXPixelsNum()); - auto det_height_pxl = static_cast(experiment.GetYPixelsNum()); + const auto geom = experiment.GetDiffractionGeometry(); + const auto det_width_pxl = static_cast(experiment.GetXPixelsNum()); + const auto det_height_pxl = static_cast(experiment.GetYPixelsNum()); - float one_over_dmax = 1.0f / settings.high_res_A; - float one_over_dmax_sq = one_over_dmax * one_over_dmax; + const float one_over_dmax = 1.0f / settings.high_res_A; + const float one_over_dmax_sq = one_over_dmax * one_over_dmax; float one_over_wavelength = 1.0f / geom.GetWavelength_A(); - Coord Astar = lattice.Astar(); - Coord Bstar = lattice.Bstar(); - Coord Cstar = lattice.Cstar(); - Coord S0 = geom.GetScatteringVector(); + const Coord Astar = lattice.Astar(); + const Coord Bstar = lattice.Bstar(); + const Coord Cstar = lattice.Cstar(); + const Coord S0 = geom.GetScatteringVector(); std::vector rot = geom.GetPoniRotMatrix().transpose().arr(); diff --git a/image_analysis/bragg_prediction/BraggPrediction.h b/image_analysis/bragg_prediction/BraggPrediction.h index e616aa0f..cf067f01 100644 --- a/image_analysis/bragg_prediction/BraggPrediction.h +++ b/image_analysis/bragg_prediction/BraggPrediction.h @@ -12,9 +12,17 @@ struct BraggPredictionSettings { float high_res_A = 1.5; - float ewald_dist_cutoff = 0.0005; int max_hkl = 100; char centering = 'P'; + + // Still parameters + float ewald_dist_cutoff = 0.0005; + + // Rotation parameters + Coord rotation_axis = Coord(1, 0, 0); + int image_number = 0; + float wedge_size_deg = 0.1f; + float mosaicity_deg = 0.1f; }; class BraggPrediction { @@ -30,11 +38,5 @@ public: const std::vector &GetReflections() const; }; -std::vector CalcBraggPredictions(const DiffractionExperiment &experiment, - const CrystalLattice &lattice, - float high_res_A = 1.5, - float ewald_dist_cutoff = 0.0005, - int max_hkl = 100, - int max_reflections = 10000); #endif //JFJOCH_BRAGGPREDICTION_H diff --git a/image_analysis/bragg_prediction/BraggPredictionRotation.cpp b/image_analysis/bragg_prediction/BraggPredictionRotation.cpp index 3c6f9bde..94a849c3 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRotation.cpp +++ b/image_analysis/bragg_prediction/BraggPredictionRotation.cpp @@ -11,6 +11,9 @@ #include "../../common/JFJochException.h" #include "../bragg_integration/SystematicAbsence.h" +BraggPredictionRotation::BraggPredictionRotation(int max_reflections) + : BraggPrediction(max_reflections) {} + namespace { inline float deg_to_rad(float deg) { return deg * (static_cast(M_PI) / 180.0f); @@ -63,12 +66,9 @@ namespace { } } // namespace -std::vector BraggPredictionRotation::Calc(const DiffractionExperiment& experiment, - const CrystalLattice& lattice, - const BraggPredictionRotationSettings& settings) const { - std::vector out; - out.reserve(200000); // tune - +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, @@ -76,134 +76,145 @@ std::vector BraggPredictionRotation::Calc(const DiffractionExperimen 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()); - // Determine prediction interval in spindle angle. - // We treat each image as [angle(image), angle(image)+wedge) - const float start_deg = gon.GetStart_deg(); - const float inc_deg = gon.GetIncrement_deg(); - const float wedge_deg = gon.GetWedge_deg(); - - float phi0_deg = gon.GetAngle_deg(static_cast(settings.image_first)); - float phi1_deg = gon.GetAngle_deg(static_cast(settings.image_last)) + wedge_deg; - - if (settings.pad_one_wedge) { - phi0_deg -= wedge_deg; - phi1_deg += wedge_deg; - } - - const float phi0 = deg_to_rad(phi0_deg); - const float phi1 = deg_to_rad(phi1_deg); - - // Resolution cutoff in reciprocal space const float one_over_dmax = 1.0f / settings.high_res_A; const float one_over_dmax_sq = one_over_dmax * one_over_dmax; - // S0 has length 1/lambda (your convention) - const Coord S0 = geom.GetScatteringVector(); - const float k2 = (S0 * S0); // (1/lambda)^2 - - // Rotation axis in lab frame (already normalized in ctor, but keep safe) - const Coord w = gon.GetAxis().Normalize(); + 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(); const Coord Cstar = lattice.Cstar(); + const Coord S0 = geom.GetScatteringVector(); - const float det_w = static_cast(experiment.GetXPixelsNum()); - const float det_h = static_cast(experiment.GetYPixelsNum()); + const std::vector rot = geom.GetPoniRotMatrix().transpose().arr(); + + // Precompute detector geometry constants + const float beam_x = geom.GetBeamX_pxl(); + 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 start_deg = gon.GetStart_deg(); + const float inc_deg = gon.GetIncrement_deg(); + 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); + + const Coord w = gon.GetAxis().Normalize(); + + int i = 0; for (int h = -settings.max_hkl; h <= settings.max_hkl; ++h) { - const Coord Ah = Astar * static_cast(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) { - const Coord AhBk = Ah + Bstar * static_cast(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) { if (systematic_absence(h, k, l, settings.centering)) continue; - const Coord g0 = AhBk + Cstar * static_cast(l); - const float g02 = g0 * g0; - - if (!(g02 > 0.0f)) - continue; - if (g02 > one_over_dmax_sq) + if (i >= max_reflections) continue; - // g0 = g_par + g_perp wrt spindle axis - const float g_par_s = g0 * w; - const Coord g_par = w * g_par_s; - const Coord g_perp = g0 - g_par; + const float g0_x = AhBk_x + Cstar.x * l; + const float g0_y = AhBk_y + Cstar.y * l; + const float g0_z = AhBk_z + Cstar.z * l; + const float g02 = g0_x * g0_x + g0_y * g0_y + g0_z * g0_z; - const float g_perp2 = g_perp * g_perp; - if (g_perp2 < 1e-12f) { - // Rotation does not move this reciprocal vector: skip for now. - // (Could be handled by checking whether it satisfies Ewald at all.) + if (g02 <= 0.0f || g02 > one_over_dmax_sq) continue; - } - // Equation: |S0 + g(phi)|^2 = |S0|^2 - // g(phi) = g_par + cos(phi) g_perp + sin(phi) (w x g_perp) - const Coord p = S0 + g_par; - const Coord w_x_gperp = w % g_perp; + const float g_par_s = g0_x * w.x + g0_y * w.y + g0_z * w.z; + const float g_par_x = w.x * g_par_s; + const float g_par_y = w.y * g_par_s; + const float g_par_z = w.z * g_par_s; - const float A = 2.0f * (p * g_perp); - const float B = 2.0f * (p * w_x_gperp); - const float D = (p * p) + g_perp2 - k2; + const float g_perp_x = g0_x - g_par_x; + const float g_perp_y = g0_y - g_par_y; + const float g_perp_z = g0_z - g_par_z; + + const float g_perp2 = g_perp_x * g_perp_x + g_perp_y * g_perp_y + g_perp_z * g_perp_z; + if (g_perp2 < 1e-12f) + continue; + + const float p_x = S0.x + g_par_x; + const float p_y = S0.y + g_par_y; + const float p_z = S0.z + g_par_z; + + const float w_x_gperp_x = w.y * g_perp_z - w.z * g_perp_y; + const float w_x_gperp_y = w.z * g_perp_x - w.x * g_perp_z; + const float w_x_gperp_z = w.x * g_perp_y - w.y * g_perp_x; + + const float A_trig = 2.0f * (p_x * g_perp_x + p_y * g_perp_y + p_z * g_perp_z); + const float B_trig = 2.0f * (p_x * w_x_gperp_x + p_y * w_x_gperp_y + p_z * w_x_gperp_z); + const float D_trig = (p_x * p_x + p_y * p_y + p_z * p_z) + g_perp2 - one_over_wavelength_sq; float sols[2]{}; - const int nsol = solve_trig(A, B, D, phi0, phi1, sols); - if (nsol == 0) - continue; + const int nsol = solve_trig(A_trig, B_trig, D_trig, phi0, phi1, 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); - // Rotate g0 by phi about w - const RotMatrix R(phi, w); - const Coord g = R * g0; - const Coord S = S0 + g; + // Rodrigues' rotation formula: g = g_par + cos(phi)*g_perp + sin(phi)*(w x g_perp) + const float g_x = g_par_x + cos_p * g_perp_x + sin_p * w_x_gperp_x; + const float g_y = g_par_y + cos_p * g_perp_y + sin_p * w_x_gperp_y; + const float g_z = g_par_z + cos_p * g_perp_z + sin_p * w_x_gperp_z; - // Project to detector using canonical geometry code - const auto [x, y] = geom.RecipToDector(g); - if (!std::isfinite(x) || !std::isfinite(y)) - continue; - if (x < 0.0f || x >= det_w || y < 0.0f || y >= det_h) + const float S_x = S0.x + g_x; + const float S_y = S0.y + g_y; + const float S_z = S0.z + g_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; + + if (S_rot_z <= 0) continue; - // Convert phi to fractional image number + 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; + + if (x < 0.0f || x >= det_width_pxl || y < 0.0f || y >= det_height_pxl) + continue; + + float d = 1.0f / sqrtf(g02); const float phi_deg = rad_to_deg(phi); - - Reflection r{}; - r.h = h; - r.k = k; - r.l = l; - r.angle_deg = phi_deg; - r.image_number = phi_deg_to_image_number(phi_deg, start_deg, inc_deg); - r.predicted_x = x; - r.predicted_y = y; - r.d = 1.0f / std::sqrt(g02); - - // diagnostic: should be ~0 - r.dist_ewald = S.Length() - std::sqrt(k2); - - r.I = 0.0f; - r.bkg = 0.0f; - r.sigma = 0.0f; - - out.push_back(r); + reflections[i] = Reflection{ + .h = h, + .k = k, + .l = l, + .angle_deg = phi_deg, + .predicted_x = x, + .predicted_y = y, + .d = d, + .dist_ewald = sqrtf(S_x * S_x + S_y * S_y + S_z * S_z) - one_over_wavelength, + }; + ++i; } } } } - - std::sort(out.begin(), out.end(), [](const Reflection &a, const Reflection &b) { - if (a.angle_deg != b.angle_deg) - return a.angle_deg < b.angle_deg; - if (a.h != b.h) return a.h < b.h; - if (a.k != b.k) return a.k < b.k; - return a.l < b.l; - }); - - return out; -} + return i; +} \ No newline at end of file diff --git a/image_analysis/bragg_prediction/BraggPredictionRotation.h b/image_analysis/bragg_prediction/BraggPredictionRotation.h index 417f6ec1..6c4de590 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRotation.h +++ b/image_analysis/bragg_prediction/BraggPredictionRotation.h @@ -10,28 +10,16 @@ #include "../../common/CrystalLattice.h" #include "../../common/DiffractionExperiment.h" #include "../../common/Reflection.h" +#include "BraggPrediction.h" -struct BraggPredictionRotationSettings { - float high_res_A = 1.5f; - int max_hkl = 100; - char centering = 'P'; - // Predict only in this image interval (inclusive) - int64_t image_first = 0; - int64_t image_last = 0; - - // If true, expands [phi_start, phi_end] by one wedge on each side - // (useful if you later want to model partials / edge effects). - bool pad_one_wedge = false; -}; - -class BraggPredictionRotation { +class BraggPredictionRotation : public BraggPrediction { public: - BraggPredictionRotation() = default; + BraggPredictionRotation(int max_reflections); - std::vector Calc(const DiffractionExperiment& experiment, - const CrystalLattice& lattice, - const BraggPredictionRotationSettings& settings) const; + int Calc(const DiffractionExperiment &experiment, + const CrystalLattice &lattice, + const BraggPredictionSettings &settings) override; };