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

BraggPredictionRotation: New implementation - looks OK, but accuracy is terrible (+/- 10 images wrong)

Browse Source
This commit is contained in:
leonarski_f
2026-01-22 14:13:52 +01:00
parent d3f311dc97
commit 45ca48cc01
4 changed files with 87 additions and 155 deletions
Download Patch File Download Diff File Expand all files Collapse all files

13
image_analysis/RotationIndexer.cpp
View File

@@ -6,6 +6,7 @@
#include "geom_refinement/XtalOptimizer.h"
#include "indexing/FFTIndexer.h"
#include "lattice_search/LatticeSearch.h"
#include <iostream>
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<RotationIndexerResult> RotationIndexer::GetLattice() {
.geom = updated_geom_
};
}
const std::vector<PredictionRotationResult> &RotationIndexer::GetPredictions() const {
return predictions;
}

4
image_analysis/RotationIndexer.h
View File

@@ -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<SpotToSave> v_;
std::vector<Coord> coords_;
@@ -51,11 +53,13 @@ class RotationIndexer {
std::optional<CrystalLattice> indexed_lattice;
std::vector<PredictionRotationResult> predictions;
void TryIndex();
public:
RotationIndexer(const DiffractionExperiment& x, IndexerThreadPool& indexer);
std::optional<RotationIndexerResult> ProcessImage(int64_t image, const std::vector<SpotToSave>& spots);
std::optional<RotationIndexerResult> GetLattice();
const std::vector<PredictionRotationResult> &GetPredictions() const;
};
#endif //JFJOCH_ROTATIONINDEXER_H

202
image_analysis/bragg_prediction/BraggPredictionRotation.cpp
View File

@@ -11,66 +11,11 @@
#include "../../common/JFJochException.h"
#include "../bragg_integration/SystematicAbsence.h"
BraggPredictionRotation::BraggPredictionRotation(int max_reflections)
: BraggPrediction(max_reflections) {}
std::vector<PredictionRotationResult> PredictRotationHKLs(const DiffractionExperiment &experiment,
const CrystalLattice &lattice,
const PredictionRotationSettings &settings) {
std::vector<PredictionRotationResult> ret;
namespace {
inline float deg_to_rad(float deg) {
return deg * (static_cast<float>(M_PI) / 180.0f);
}
inline float rad_to_deg(float rad) {
return rad * (180.0f / static_cast<float>(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<float>(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<float>(experiment.GetXPixelsNum());
const float det_height_pxl = static_cast<float>(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<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);
// 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;
}
std::ranges::sort(ret,
[](const PredictionRotationResult& a, const PredictionRotationResult& b) {
return a.angle_deg < b.angle_deg;
});
return ret;
}

23
image_analysis/bragg_prediction/BraggPredictionRotation.h
View File

@@ -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<PredictionRotationResult> PredictRotationHKLs(const DiffractionExperiment &experiment,
const CrystalLattice &lattice,
const PredictionRotationSettings &settings);
#endif //JFJOCH_BRAGGPREDICTIONROTATION_H