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

XtalOptimizer: Add option to refine rot1/rot2 + add bounds (+/- 3.0 deg for rot1/rot2 and 10% for detector distance)

Browse Source
This commit is contained in:
leonarski_f
2025-12-01 15:07:15 +01:00
parent 739215c5ac
commit d86a32ad43
2 changed files with 40 additions and 22 deletions
Download Patch File Download Diff File Expand all files Collapse all files

61
image_analysis/geom_refinement/XtalOptimizer.cpp
View File

@@ -9,7 +9,6 @@ struct XtalResidual {
XtalResidual(double x, double y,
double lambda,
double pixel_size,
double rot1, double rot2,
const Eigen::Matrix3d &gonio_back_rot,
double exp_h, double exp_k,
double exp_l,
@@ -20,8 +19,6 @@ struct XtalResidual {
exp_h(exp_h),
exp_k(exp_k),
exp_l(exp_l),
rot1(rot1),
rot2(rot2),
gonio_back_rot_(gonio_back_rot),
symmetry(symmetry) {
}
@@ -30,14 +27,15 @@ struct XtalResidual {
bool operator()(const T *const beam_x,
const T *const beam_y,
const T *const distance_mm,
const T *const detector_rot,
const T *const p0,
const T *const p1,
const T *const p2,
T *residual) const {
T c1 = ceres::cos(T(rot1));
T c2 = ceres::cos(T(rot2));
T s1 = ceres::sin(T(rot1));
T s2 = ceres::sin(T(rot2));
T c1 = ceres::cos(T(detector_rot[0]));
T c2 = ceres::cos(T(detector_rot[1]));
T s1 = ceres::sin(T(detector_rot[0]));
T s2 = ceres::sin(T(detector_rot[1]));
// x_lab in mm
T x_lab = (T(obs_x) - beam_x[0]) * T(pixel_size);
@@ -141,7 +139,6 @@ struct XtalResidual {
const double exp_h;
const double exp_k;
const double exp_l;
const double rot1, rot2;
const Eigen::Matrix3d gonio_back_rot_;
gemmi::CrystalSystem symmetry;
};
@@ -379,6 +376,8 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
double beam_y = data.geom.GetBeamY_pxl();
double distance_mm = data.geom.GetDetectorDistance_mm();
double detector_rot[2] = {data.geom.GetPoniRot1_rad(), data.geom.GetPoniRot2_rad()};
ceres::Problem problem;
double latt_vec0[3], latt_vec1[3], latt_vec2[3];
@@ -447,16 +446,11 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
continue;
problem.AddResidualBlock(
new ceres::AutoDiffCostFunction<XtalResidual, 3, 1, 1, 1, 3, 3, 3>(
new ceres::AutoDiffCostFunction<XtalResidual, 3, 1, 1, 1, 2, 3, 3, 3>(
new XtalResidual(pt.x, pt.y,
data.geom.GetWavelength_A(),
data.geom.GetPixelSize_mm(),
data.geom.GetPoniRot1_rad(),
data.geom.GetPoniRot2_rad(),
gonio_back_rot,
h, k, l,
data.crystal_system)),
@@ -464,6 +458,7 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
&beam_x,
&beam_y,
&distance_mm,
detector_rot,
latt_vec0,
latt_vec1,
latt_vec2
@@ -475,12 +470,27 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
if (!data.refine_distance_mm)
problem.SetParameterBlockConstant(&distance_mm);
else {
const double dist_range = 0.1;
problem.SetParameterLowerBound(&distance_mm, 0, distance_mm * (1.0 - dist_range));
problem.SetParameterLowerBound(&distance_mm, 0, distance_mm * (1.0 + dist_range));
}
if (!data.refine_beam_center) {
problem.SetParameterBlockConstant(&beam_x);
problem.SetParameterBlockConstant(&beam_y);
}
if (!data.refine_detector_angles) {
problem.SetParameterBlockConstant(detector_rot);
} else {
const double rot_range = 3.0 / 180.0 * M_PI;
for (int i = 0; i < 2; ++i) {
problem.SetParameterLowerBound(detector_rot, i, detector_rot[i] - rot_range);
problem.SetParameterUpperBound(detector_rot, i, detector_rot[i] + rot_range);
}
}
// Parameter bounds
// Lengths
for (int i = 0; i < 3; ++i) {
@@ -522,6 +532,9 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
if (data.refine_distance_mm)
data.geom.DetectorDistance_mm(distance_mm);
if (data.refine_detector_angles)
data.geom.PoniRot1_rad(detector_rot[0]).PoniRot2_rad(detector_rot[1]);
if (data.crystal_system == gemmi::CrystalSystem::Orthorhombic)
data.latt = AngleAxisAndLengthsToLattice(latt_vec0, latt_vec1, false);
else if (data.crystal_system == gemmi::CrystalSystem::Tetragonal) {
@@ -551,18 +564,22 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
const double sg = std::sin(latt_vec2[2]);
// a along x, b in x-y, c general
B(0,0) = 1.0; B(1,0) = 0.0; B(2,0) = 0.0;
B(0,1) = cg; B(1,1) = sg; B(2,1) = 0.0;
B(0, 0) = 1.0;
B(1, 0) = 0.0;
B(2, 0) = 0.0;
B(0, 1) = cg;
B(1, 1) = sg;
B(2, 1) = 0.0;
const double cx = cb;
const double cy = (ca - cb*cg) / sg;
const double cz = std::sqrt(std::max(0.0, 1.0 - cx*cx - cy*cy));
const double cy = (ca - cb * cg) / sg;
const double cz = std::sqrt(std::max(0.0, 1.0 - cx * cx - cy * cy));
B(0,2) = cx;
B(1,2) = cy;
B(2,2) = cz;
B(0, 2) = cx;
B(1, 2) = cy;
B(2, 2) = cz;
Eigen::DiagonalMatrix<double,3> D(latt_vec1[0], latt_vec1[1], latt_vec1[2]);
Eigen::DiagonalMatrix<double, 3> D(latt_vec1[0], latt_vec1[1], latt_vec1[2]);
Eigen::Matrix3d latt = R * D * B;
data.latt = CrystalLattice(Coord(latt(0, 0), latt(1, 0), latt(2, 0)),

1
image_analysis/geom_refinement/XtalOptimizer.h
View File

@@ -25,6 +25,7 @@ struct XtalOptimizerData {
bool refine_beam_center = true;
bool refine_distance_mm = false;
bool refine_detector_angles = false;
std::optional<GoniometerAxis> axis;