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

IndexAndRefine: Use small angle approximation for optimizing rotation data

Browse Source
This commit is contained in:
leonarski_f
2026-01-30 14:44:54 +01:00
parent 52339a7f85
commit a5e582574e
4 changed files with 84 additions and 49 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/geom_refinement/SimpleRotXtalOptimizer.cpp
+42 -14
View File
@@ -8,7 +8,8 @@
static bool SimpleRotXtalOptimizerInternal(
SimpleRotXtalOptimizerData &data,
const std::vector<SpotToSave> &spots,
float tolerance)
float tolerance,
Eigen::Matrix3d &r_total)
{
Coord vec0 = data.latt.Vec0();
Coord vec1 = data.latt.Vec1();
@@ -59,10 +60,11 @@ static bool SimpleRotXtalOptimizerInternal(
const auto &h = obs[i];
const auto &e = exp[i];
// Cross-product matrix of h_obs
A.row(3*i + 0) << 0.0, -h.z(), h.y();
A.row(3*i + 1) << h.z(), 0.0, -h.x();
A.row(3*i + 2) << -h.y(), h.x(), 0.0;
// Cross-product matrix of h_obs: [h]× such that [h]× · ω = h × ω
// But we need ω × h = -h × ω, so use -[h]×
A.row(3*i + 0) << 0.0, h.z(), -h.y();
A.row(3*i + 1) << -h.z(), 0.0, h.x();
A.row(3*i + 2) << h.y(), -h.x(), 0.0;
b.segment<3>(3*i) = e - h;
}
@@ -74,21 +76,47 @@ static bool SimpleRotXtalOptimizerInternal(
if (!omega.allFinite())
return false;
if (omega.norm() > 0.05) {
return false;
// > ~3 degrees — warn or reject
}
// Apply incremental rotation and accumulate it
const double angle = omega.norm();
Eigen::Matrix3d r_inc = Eigen::Matrix3d::Identity();
if (angle > 0.0)
r_inc = Eigen::AngleAxisd(angle, omega / angle).toRotationMatrix();
// Apply small-angle rotation to lattice
// TODO: data.latt.RotateSmallAngle(omega);
r_total = r_inc * r_total;
if (angle > 0.0) {
// Build the current lattice matrix (columns = a, b, c)
Eigen::Matrix3d L;
L.col(0) = Eigen::Vector3d(vec0.x, vec0.y, vec0.z);
L.col(1) = Eigen::Vector3d(vec1.x, vec1.y, vec1.z);
L.col(2) = Eigen::Vector3d(vec2.x, vec2.y, vec2.z);
// The HKL-space rotation R transforms: h_new = R * h_old
// This corresponds to: L_new = L_old * R^T (real-space lattice)
Eigen::Matrix3d L_new = L * r_inc.transpose();
data.latt = CrystalLattice(
Coord(L_new(0,0), L_new(1,0), L_new(2,0)),
Coord(L_new(0,1), L_new(1,1), L_new(2,1)),
Coord(L_new(0,2), L_new(1,2), L_new(2,2))
);
}
return true;
}
bool SimpleRotXtalOptimizer(SimpleRotXtalOptimizerData &data, const std::vector<SpotToSave> &spots) {
Eigen::Matrix3d r_total = Eigen::Matrix3d::Identity();
if (!SimpleRotXtalOptimizerInternal(data, spots, 0.3))
if (!SimpleRotXtalOptimizerInternal(data, spots, 0.3, r_total))
return false;
SimpleRotXtalOptimizerInternal(data, spots, 0.2);
return SimpleRotXtalOptimizerInternal(data, spots, 0.1);
SimpleRotXtalOptimizerInternal(data, spots, 0.2, r_total);
const bool ok = SimpleRotXtalOptimizerInternal(data, spots, 0.1, r_total);
Eigen::AngleAxisd aa(r_total);
data.rotation[0] = aa.axis().x();
data.rotation[1] = aa.axis().y();
data.rotation[2] = aa.axis().z();
data.angle = aa.angle() * 180.0 / M_PI;
return ok;
}