From 9e01f73cd489353d1ff2b9b017fb01c2c2b0f9f7 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Tue, 20 Jan 2026 13:07:14 +0100 Subject: [PATCH] BraggPredictionRotation: Work in progress --- .../BraggPredictionRotation.cpp | 66 +++++++++---------- 1 file changed, 33 insertions(+), 33 deletions(-) diff --git a/image_analysis/bragg_prediction/BraggPredictionRotation.cpp b/image_analysis/bragg_prediction/BraggPredictionRotation.cpp index 40b90114..5c55ded4 100644 --- a/image_analysis/bragg_prediction/BraggPredictionRotation.cpp +++ b/image_analysis/bragg_prediction/BraggPredictionRotation.cpp @@ -103,8 +103,6 @@ int BraggPredictionRotation::Calc(const DiffractionExperiment &experiment, 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)); @@ -117,7 +115,7 @@ int BraggPredictionRotation::Calc(const DiffractionExperiment &experiment, const float mosaicity_rad = deg_to_rad(settings.mosaicity_deg); const float half_mosaicity = mosaicity_rad * 0.5f; - const Coord w = gon.GetAxis().Normalize(); + const Coord m2 = gon.GetAxis().Normalize(); int i = 0; @@ -138,38 +136,40 @@ int BraggPredictionRotation::Calc(const DiffractionExperiment &experiment, if (i >= max_reflections) continue; - const float recip_x = AhBk_x + Cstar.x * l; - const float recip_y = AhBk_y + Cstar.y * l; - const float recip_z = AhBk_z + Cstar.z * l; - const float recip_sq = recip_x * recip_x + recip_y * recip_y + recip_z * recip_z; + 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; - if (recip_sq <= 0.0f || recip_sq > one_over_dmax_sq) + if (p0_sq <= 0.0f || p0_sq > one_over_dmax_sq) continue; - const float g_par_s = recip_x * w.x + recip_y * w.y + recip_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 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 g_perp_x = recip_x - g_par_x; - const float g_perp_y = recip_y - g_par_y; - const float g_perp_z = recip_z - g_par_z; + 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 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) + 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) 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 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; - 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 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; - 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; + // 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); @@ -183,13 +183,13 @@ int BraggPredictionRotation::Calc(const DiffractionExperiment &experiment, const float sin_p = sinf(phi); // 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; + 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 + g_x; - const float S_y = S0.y + g_y; - const float S_z = S0.z + g_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; @@ -206,7 +206,7 @@ int BraggPredictionRotation::Calc(const DiffractionExperiment &experiment, 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(recip_sq); + float d = 1.0f / sqrtf(p0_sq); const float phi_deg = rad_to_deg(phi); reflections[i] = Reflection{ .h = h,