mx/Jungfraujoch
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 42 Wiki Activity

BraggPredictionRotation: Work in progress

Browse Source
This commit is contained in:
leonarski_f
2026-01-20 13:07:14 +01:00
parent 2e29834fb8
commit 9e01f73cd4
1 changed files with 33 additions and 33 deletions
Download Patch File Download Diff File Expand all files Collapse all files

66
image_analysis/bragg_prediction/BraggPredictionRotation.cpp
View File

@@ -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<float>(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,