BraggPrediction: Work in progress

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);

image_analysis/bragg_prediction/BraggPrediction.cpp

@@ -14,19 +14,19 @@ const std::vector<Reflection> &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<float>(experiment.GetXPixelsNum());
-    auto det_height_pxl = static_cast<float>(experiment.GetYPixelsNum());
+    const auto geom = experiment.GetDiffractionGeometry();
+    const auto det_width_pxl = static_cast<float>(experiment.GetXPixelsNum());
+    const auto det_height_pxl = static_cast<float>(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<float> rot = geom.GetPoniRotMatrix().transpose().arr();
 

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<Reflection> &GetReflections() const;
 };
 
-std::vector<Reflection> 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

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<float>(M_PI) / 180.0f);
@@ -63,12 +66,9 @@ namespace {
     }
 } // namespace
 
-std::vector<Reflection> BraggPredictionRotation::Calc(const DiffractionExperiment& experiment,
-                                                        const CrystalLattice& lattice,
-                                                        const BraggPredictionRotationSettings& settings) const {
-    std::vector<Reflection> 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<Reflection> BraggPredictionRotation::Calc(const DiffractionExperimen
 
     const GoniometerAxis& gon = *gon_opt;
     const auto geom = experiment.GetDiffractionGeometry();
+    const float det_width_pxl = static_cast<float>(experiment.GetXPixelsNum());
+    const float det_height_pxl = static_cast<float>(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<float>(settings.image_first));
-    float phi1_deg = gon.GetAngle_deg(static_cast<float>(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<float>(experiment.GetXPixelsNum());
-    const float det_h = static_cast<float>(experiment.GetYPixelsNum());
+    const std::vector<float> 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<float>(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<float>(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<float>(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<float>(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;
+}

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<Reflection> Calc(const DiffractionExperiment& experiment,
-                                 const CrystalLattice& lattice,
-                                 const BraggPredictionRotationSettings& settings) const;
+    int Calc(const DiffractionExperiment &experiment,
+             const CrystalLattice &lattice,
+             const BraggPredictionSettings &settings) override;
 };