leonarski_fandClaude Opus 4.8
acd2025676
The azimuthal-integration solid-angle correction used cos^3(2*theta), where 2*theta is the true scattering angle (from LabCoord, including detector tilt). The solid angle of a flat pixel actually depends on the incidence angle to the detector normal, cos(alpha) = det_distance / |detector-frame position|, which is invariant under detector tilt (rot1/rot2/rot3). Only for an untilted detector do the two agree. Switch CalcAzIntSolidAngleCorr(x,y) to the tilt-invariant form, matching PyFAI solidAngleArray and MAX IV azint. Drop the q-only overload (it can only ever be the untilted approximation and was used only in tests) and move its test onto the (x,y) form; add a tilt-invariance test. XtalOptimizer's residual reconstructed each spot's lab position from rot1/rot2 only, hardcoding rot3 = 0, while the rest of the pipeline (and its own spot selection) used the full PONI rotation. An imported non-zero rot3 was therefore silently dropped during refinement. Bake rot3 into the residual as a fixed Rz(-rot3) so refinement stays consistent (no-op when rot3 == 0). Polarization and azimuthal binning already honoured rot3 via the full PONI rotation (Phi_rad), validated against PyFAI chi() by the existing rot3 phi tests. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
627 lines
25 KiB
C++
627 lines
25 KiB
C++
// SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
|
||
// SPDX-License-Identifier: GPL-3.0-only
|
||
|
||
#include "../../common/JFJochMath.h"
|
||
#include <Eigen/Dense>
|
||
|
||
#include "XtalOptimizer.h"
|
||
#include "ceres/ceres.h"
|
||
#include "ceres/rotation.h"
|
||
#include "LatticeReduction.h"
|
||
|
||
struct XtalResidual {
|
||
XtalResidual(double x, double y,
|
||
double lambda,
|
||
double pixel_size,
|
||
double rot3,
|
||
double angle_rad,
|
||
double exp_h, double exp_k,
|
||
double exp_l,
|
||
gemmi::CrystalSystem symmetry)
|
||
: obs_x(x), obs_y(y),
|
||
inv_lambda(1.0/lambda),
|
||
pixel_size(pixel_size),
|
||
rot3(rot3),
|
||
exp_h(exp_h),
|
||
exp_k(exp_k),
|
||
exp_l(exp_l),
|
||
angle_rad(angle_rad),
|
||
symmetry(symmetry) {
|
||
if (std::fabs(lambda) < 1e-6)
|
||
throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
|
||
"Lambda cannot be close to zero");
|
||
}
|
||
|
||
template<typename T>
|
||
bool operator()(const T *const beam,
|
||
const T *const distance_mm,
|
||
const T *const detector_rot,
|
||
const T *const rotation_axis,
|
||
const T *const p0,
|
||
const T *const p1,
|
||
const T *const p2,
|
||
T *residual) const {
|
||
// PyFAI convention: poni_rot = Rz(-rot3) * Rx(-rot2) * Ry(+rot1).
|
||
// detector_rot[0] = rot1, detector_rot[1] = rot2 are refined; rot3 is fixed
|
||
// (e.g. from a PONI import) and baked in here as a constant so that a non-zero
|
||
// rot3 is not silently dropped during refinement.
|
||
const T rot1 = detector_rot[0];
|
||
const T rot2 = detector_rot[1];
|
||
|
||
// Ry(+rot1): rotation around Y-axis
|
||
const T c1 = ceres::cos(rot1);
|
||
const T s1 = ceres::sin(rot1);
|
||
|
||
// Rx(-rot2): rotation around X-axis with inverted sign (PyFAI left-handed)
|
||
const T c2 = ceres::cos(rot2);
|
||
const T s2 = ceres::sin(rot2);
|
||
|
||
// Rz(-rot3): rotation around Z (beam); constant, identity when rot3 == 0
|
||
const T c3 = T(cos(rot3));
|
||
const T s3 = T(sin(rot3));
|
||
|
||
// Detector coordinates in mm
|
||
const T det_x = (T(obs_x) - beam[0]) * T(pixel_size);
|
||
const T det_y = (T(obs_y) - beam[1]) * T(pixel_size);
|
||
const T det_z = T(distance_mm[0]);
|
||
|
||
// Apply Ry(rot1) first: rotate around Y
|
||
const T t1_x = c1 * det_x + s1 * det_z;
|
||
const T t1_y = det_y;
|
||
const T t1_z = -s1 * det_x + c1 * det_z;
|
||
|
||
// Then apply Rx(-rot2): rotate around X
|
||
const T t2_x = t1_x;
|
||
const T t2_y = c2 * t1_y + s2 * t1_z;
|
||
const T t2_z = -s2 * t1_y + c2 * t1_z;
|
||
|
||
// Then apply Rz(-rot3): rotate around Z (beam)
|
||
const T x = c3 * t2_x + s3 * t2_y;
|
||
const T y = -s3 * t2_x + c3 * t2_y;
|
||
const T z = t2_z;
|
||
|
||
// convert to recip space
|
||
const T lab_norm = ceres::sqrt(x * x + y * y + z * z);
|
||
const T inv_norm = T(1) / lab_norm;
|
||
|
||
T recip_raw[3];
|
||
recip_raw[0] = x * inv_norm * T(inv_lambda);
|
||
recip_raw[1] = y * inv_norm * T(inv_lambda);
|
||
recip_raw[2] = (z * inv_norm - T(1.0)) * T(inv_lambda);
|
||
|
||
// Apply goniometer "back-to-start" rotation:
|
||
// brings observed reciprocal from image orientation into reference crystal frame
|
||
const T aa_back[3] = {
|
||
T(angle_rad) * rotation_axis[0],
|
||
T(angle_rad) * rotation_axis[1],
|
||
T(angle_rad) * rotation_axis[2]
|
||
};
|
||
|
||
T recip_obs[3];
|
||
ceres::AngleAxisRotatePoint(aa_back, recip_raw, recip_obs);
|
||
|
||
const Eigen::Matrix<T, 3, 1> e_obs_recip(recip_obs[0], recip_obs[1], recip_obs[2]);
|
||
|
||
// Build unit cell lengths and B (convention: columns are a, b, c prior to global rotation)
|
||
Eigen::Matrix<T, 3, 1> e_uc_len = Eigen::Matrix<T, 3, 1>::Zero();
|
||
Eigen::Matrix<T, 3, 3> B = Eigen::Matrix<T, 3, 3>::Identity();
|
||
|
||
if (symmetry == gemmi::CrystalSystem::Hexagonal) {
|
||
e_uc_len << p1[0], p1[0], p1[2];
|
||
B(0, 1) = T(-0.5); // cos(120)
|
||
B(1, 1) = T(sqrt(3.0) / 2.0); // sin(120)
|
||
} else if (symmetry == gemmi::CrystalSystem::Orthorhombic) {
|
||
e_uc_len << p1[0], p1[1], p1[2];
|
||
} else if (symmetry == gemmi::CrystalSystem::Tetragonal) {
|
||
e_uc_len << p1[0], p1[0], p1[2];
|
||
} else if (symmetry == gemmi::CrystalSystem::Cubic) {
|
||
e_uc_len << p1[0], p1[0], p1[0];
|
||
} else if (symmetry == gemmi::CrystalSystem::Monoclinic) {
|
||
// Unique axis b: alpha = gamma = 90°, beta free (angle between a and c)
|
||
e_uc_len << p1[0], p1[1], p1[2];
|
||
B(0, 2) = ceres::cos(p2[0]);
|
||
B(2, 2) = ceres::sin(p2[0]);
|
||
} else {
|
||
// Triclinic: p1 = (a,b,c), p2 = (alpha, beta, gamma) in radians
|
||
const T ca = ceres::cos(p2[0]);
|
||
const T cb = ceres::cos(p2[1]);
|
||
const T cg = ceres::cos(p2[2]);
|
||
const T sg = ceres::sin(p2[2]);
|
||
|
||
e_uc_len << p1[0], p1[1], p1[2];
|
||
|
||
B(0, 0) = T(1);
|
||
B(1, 0) = T(0);
|
||
B(2, 0) = T(0);
|
||
B(0, 1) = cg;
|
||
B(1, 1) = sg;
|
||
B(2, 1) = T(0);
|
||
|
||
// c vector components:
|
||
const T cx = cb;
|
||
const T cy = (ca - cb * cg) / sg;
|
||
const T v = T(1) - cx * cx - cy * cy;
|
||
const T cz = (v >= T(0)) ? ceres::sqrt(v) : T(0);
|
||
|
||
B(0, 2) = cx;
|
||
B(1, 2) = cy;
|
||
B(2, 2) = cz;
|
||
}
|
||
|
||
// Build unrotated direct lattice columns: (B * D), then rotate them by p0.
|
||
// This avoids AngleAxisToRotationMatrix + matrix multiplications.
|
||
const T L0 = e_uc_len[0];
|
||
const T L1 = e_uc_len[1];
|
||
const T L2 = e_uc_len[2];
|
||
|
||
T col0_unrot[3] = {B(0, 0) * L0, B(1, 0) * L0, B(2, 0) * L0};
|
||
T col1_unrot[3] = {B(0, 1) * L1, B(1, 1) * L1, B(2, 1) * L1};
|
||
T col2_unrot[3] = {B(0, 2) * L2, B(1, 2) * L2, B(2, 2) * L2};
|
||
|
||
T col0_rot[3], col1_rot[3], col2_rot[3];
|
||
ceres::AngleAxisRotatePoint(p0, col0_unrot, col0_rot);
|
||
ceres::AngleAxisRotatePoint(p0, col1_unrot, col1_rot);
|
||
ceres::AngleAxisRotatePoint(p0, col2_unrot, col2_rot);
|
||
|
||
const Eigen::Matrix<T, 3, 1> A(col0_rot[0], col0_rot[1], col0_rot[2]);
|
||
const Eigen::Matrix<T, 3, 1> Bv(col1_rot[0], col1_rot[1], col1_rot[2]);
|
||
const Eigen::Matrix<T, 3, 1> C(col2_rot[0], col2_rot[1], col2_rot[2]);
|
||
|
||
const Eigen::Matrix<T, 3, 1> BxC = Bv.cross(C);
|
||
const Eigen::Matrix<T, 3, 1> CxA = C.cross(A);
|
||
const Eigen::Matrix<T, 3, 1> AxB = A.cross(Bv);
|
||
|
||
const T V = A.dot(BxC);
|
||
const T invV = T(1) / V;
|
||
|
||
const Eigen::Matrix<T, 3, 1> Astar = BxC * invV;
|
||
const Eigen::Matrix<T, 3, 1> Bstar = CxA * invV;
|
||
const Eigen::Matrix<T, 3, 1> Cstar = AxB * invV;
|
||
|
||
const T h = T(exp_h);
|
||
const T k = T(exp_k);
|
||
const T l = T(exp_l);
|
||
|
||
const Eigen::Matrix<T, 3, 1> e_pred_recip = Astar * h + Bstar * k + Cstar * l;
|
||
|
||
residual[0] = e_obs_recip[0] - e_pred_recip[0];
|
||
residual[1] = e_obs_recip[1] - e_pred_recip[1];
|
||
residual[2] = e_obs_recip[2] - e_pred_recip[2];
|
||
|
||
return true;
|
||
}
|
||
|
||
const double obs_x, obs_y;
|
||
const double inv_lambda;
|
||
const double pixel_size;
|
||
const double rot3;
|
||
const double exp_h;
|
||
const double exp_k;
|
||
const double exp_l;
|
||
const double angle_rad;
|
||
gemmi::CrystalSystem symmetry;
|
||
};
|
||
|
||
struct XtalResidualRotationOnlyPrecomp {
|
||
XtalResidualRotationOnlyPrecomp(const Coord &recip_obs,
|
||
const CrystalLattice &latt,
|
||
double h, double k, double l)
|
||
: s_obs(recip_obs),
|
||
astar(latt.Astar()), bstar(latt.Bstar()), cstar(latt.Cstar()),
|
||
h(h), k(k), l(l) {
|
||
}
|
||
|
||
template<typename T>
|
||
bool operator()(const T *const rot_aa, T *residual) const {
|
||
const T astar_unrot[3] = {T(astar.x), T(astar.y), T(astar.z)};
|
||
const T bstar_unrot[3] = {T(bstar.x), T(bstar.y), T(bstar.z)};
|
||
const T cstar_unrot[3] = {T(cstar.x), T(cstar.y), T(cstar.z)};
|
||
|
||
T astar_rot[3], bstar_rot[3], cstar_rot[3];
|
||
|
||
ceres::AngleAxisRotatePoint(rot_aa, astar_unrot, astar_rot);
|
||
ceres::AngleAxisRotatePoint(rot_aa, bstar_unrot, bstar_rot);
|
||
ceres::AngleAxisRotatePoint(rot_aa, cstar_unrot, cstar_rot);
|
||
|
||
const Eigen::Matrix<T, 3, 1> s_pred(T(h) * astar_rot[0] + T(k) * bstar_rot[0] + T(l) * cstar_rot[0],
|
||
T(h) * astar_rot[1] + T(k) * bstar_rot[1] + T(l) * cstar_rot[1],
|
||
T(h) * astar_rot[2] + T(k) * bstar_rot[2] + T(l) * cstar_rot[2]
|
||
);
|
||
|
||
// Residual in reciprocal space
|
||
residual[0] = T(s_obs.x) - s_pred[0];
|
||
residual[1] = T(s_obs.y) - s_pred[1];
|
||
residual[2] = T(s_obs.z) - s_pred[2];
|
||
return true;
|
||
}
|
||
|
||
const Coord s_obs;
|
||
const Coord astar, bstar, cstar;
|
||
const double h, k, l;
|
||
};
|
||
|
||
// Regularizer: penalises ||rot_aa|| to prefer the smallest rotation that
|
||
// explains the data. Weight should be chosen in the same units as the
|
||
// reciprocal-space residuals (Å⁻¹ per radian). A value of ~0.01–0.1 is
|
||
// typically enough to break degeneracy without biasing the solution.
|
||
struct RotationNormRegularizer {
|
||
explicit RotationNormRegularizer(double weight) : weight(weight) {}
|
||
|
||
template<typename T>
|
||
bool operator()(const T *const rot_aa, T *residual) const {
|
||
residual[0] = T(weight) * rot_aa[0];
|
||
residual[1] = T(weight) * rot_aa[1];
|
||
residual[2] = T(weight) * rot_aa[2];
|
||
return true;
|
||
}
|
||
|
||
const double weight;
|
||
};
|
||
|
||
bool XtalOptimizerInternal(XtalOptimizerData &data,
|
||
const std::vector<std::vector<SpotToSave>> &spots,
|
||
const float tolerance) {
|
||
try {
|
||
Coord vec0 = data.latt.Vec0();
|
||
Coord vec1 = data.latt.Vec1();
|
||
Coord vec2 = data.latt.Vec2();
|
||
double beta = data.latt.GetUnitCell().beta;
|
||
|
||
// Initial guess for the parameters
|
||
double beam[2] = {data.geom.GetBeamX_pxl(), 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] = {0.0, 0.0, 0.0};
|
||
double latt_vec1[3] = {0.0, 0.0, 0.0};
|
||
double latt_vec2[3] = {0.0, 0.0, 0.0};
|
||
|
||
double rot_vec[3] = {1, 0, 0};
|
||
|
||
switch (data.crystal_system) {
|
||
case gemmi::CrystalSystem::Orthorhombic:
|
||
LatticeToRodriguesAndLengths_GS(data.latt, latt_vec0, latt_vec1);
|
||
break;
|
||
case gemmi::CrystalSystem::Tetragonal:
|
||
LatticeToRodriguesAndLengths_GS(data.latt, latt_vec0, latt_vec1);
|
||
latt_vec1[0] = (latt_vec1[0] + latt_vec1[1]) / 2.0;
|
||
break;
|
||
case gemmi::CrystalSystem::Cubic:
|
||
LatticeToRodriguesAndLengths_GS(data.latt, latt_vec0, latt_vec1);
|
||
latt_vec1[0] = (latt_vec1[0] + latt_vec1[1] + latt_vec1[2]) / 3.0;
|
||
break;
|
||
case gemmi::CrystalSystem::Hexagonal:
|
||
LatticeToRodriguesAndLengths_Hex(data.latt, latt_vec0, latt_vec1);
|
||
break;
|
||
case gemmi::CrystalSystem::Monoclinic:
|
||
LatticeToRodriguesLengthsBeta_Mono(data.latt, latt_vec0, latt_vec1, beta);
|
||
latt_vec2[0] = beta;
|
||
latt_vec2[1] = 0.0;
|
||
latt_vec2[2] = 0.0;
|
||
break;
|
||
default:
|
||
// Triclinic: initialize a,b,c and α,β,γ from current unit cell
|
||
LatticeToRodriguesAndLengths_GS(data.latt, latt_vec0, latt_vec1);
|
||
auto uc = data.latt.GetUnitCell();
|
||
latt_vec2[0] = uc.alpha * PI / 180.0;
|
||
latt_vec2[1] = uc.beta * PI / 180.0;
|
||
latt_vec2[2] = uc.gamma * PI / 180.0;
|
||
break;
|
||
}
|
||
|
||
if (data.axis) {
|
||
rot_vec[0] = data.axis->GetAxis().x;
|
||
rot_vec[1] = data.axis->GetAxis().y;
|
||
rot_vec[2] = data.axis->GetAxis().z;
|
||
}
|
||
|
||
const float tolerance_sq = tolerance * tolerance;
|
||
|
||
for (int i = 0; i < spots.size(); i++) {
|
||
if (spots[i].empty())
|
||
continue;
|
||
|
||
double angle_rad = 0.0;
|
||
std::optional<RotMatrix> rot_matr;
|
||
|
||
if (data.axis) {
|
||
const float angle_deg = data.axis->GetAngle_deg(i) + data.axis->GetWedge_deg() / 2.0;
|
||
angle_rad = angle_deg * PI / 180.0;
|
||
rot_matr = data.axis->GetTransformationAngle(angle_deg);
|
||
}
|
||
|
||
// Add residuals for each point
|
||
for (const auto &pt: spots[i]) {
|
||
if (!data.index_ice_rings && pt.ice_ring)
|
||
continue;
|
||
|
||
Coord recip = pt.ReciprocalCoord(data.geom);
|
||
|
||
if (rot_matr)
|
||
recip = rot_matr.value() * recip;
|
||
|
||
double h_fp = recip * vec0;
|
||
double k_fp = recip * vec1;
|
||
double l_fp = recip * vec2;
|
||
|
||
double h = std::round(h_fp);
|
||
double k = std::round(k_fp);
|
||
double l = std::round(l_fp);
|
||
|
||
double norm_sq = (h - h_fp) * (h - h_fp) + (k - k_fp) * (k - k_fp) + (l - l_fp) * (l - l_fp);
|
||
|
||
if (norm_sq > tolerance_sq)
|
||
continue;
|
||
|
||
problem.AddResidualBlock(
|
||
new ceres::AutoDiffCostFunction<XtalResidual, 3, 2, 1, 2, 3, 3, 3, 3>(
|
||
new XtalResidual(pt.x, pt.y,
|
||
data.geom.GetWavelength_A(),
|
||
data.geom.GetPixelSize_mm(),
|
||
data.geom.GetPoniRot3_rad(),
|
||
angle_rad,
|
||
h, k, l,
|
||
data.crystal_system)),
|
||
nullptr,
|
||
beam,
|
||
&distance_mm,
|
||
detector_rot,
|
||
rot_vec,
|
||
latt_vec0,
|
||
latt_vec1,
|
||
latt_vec2
|
||
);
|
||
}
|
||
}
|
||
|
||
if (problem.NumResidualBlocks() < data.min_spots)
|
||
return false;
|
||
|
||
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.SetParameterUpperBound(&distance_mm, 0, distance_mm * (1.0 + dist_range));
|
||
}
|
||
|
||
if (!data.refine_beam_center)
|
||
problem.SetParameterBlockConstant(beam);
|
||
|
||
if (!data.refine_detector_angles) {
|
||
problem.SetParameterBlockConstant(detector_rot);
|
||
} else {
|
||
const double rot_range = 3.0 / 180.0 * 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);
|
||
}
|
||
}
|
||
|
||
if (!data.refine_rotation_axis) {
|
||
problem.SetParameterBlockConstant(rot_vec);
|
||
}
|
||
|
||
if (!data.refine_unit_cell) {
|
||
problem.SetParameterBlockConstant(latt_vec1);
|
||
problem.SetParameterBlockConstant(latt_vec2);
|
||
} else {
|
||
// Parameter bounds
|
||
// Lengths
|
||
for (int i = 0; i < 3; ++i) {
|
||
problem.SetParameterLowerBound(latt_vec1, i, data.min_length_A);
|
||
problem.SetParameterUpperBound(latt_vec1, i, data.max_length_A);
|
||
}
|
||
|
||
if (data.crystal_system == gemmi::CrystalSystem::Monoclinic) {
|
||
const double beta_lo = std::max(1e-6, PI * (data.min_angle_deg / 180.0));
|
||
const double beta_hi = std::min(PI - 1e-6, PI * (data.max_angle_deg / 180.0));
|
||
problem.SetParameterLowerBound(latt_vec2, 0, beta_lo);
|
||
problem.SetParameterUpperBound(latt_vec2, 0, beta_hi);
|
||
} else if (data.crystal_system == gemmi::CrystalSystem::Triclinic) {
|
||
// α, β, γ bounds (radians)
|
||
const double alo = PI * (data.min_angle_deg / 180.0);
|
||
const double ahi = PI * (data.max_angle_deg / 180.0);
|
||
for (int i = 0; i < 3; ++i) {
|
||
problem.SetParameterLowerBound(latt_vec2, i, alo);
|
||
problem.SetParameterUpperBound(latt_vec2, i, ahi);
|
||
}
|
||
} else {
|
||
// Orthorhombic / Tetragonal / Cubic / Hexagonal:
|
||
// latt_vec2 has no meaning for these systems — always freeze it.
|
||
problem.SetParameterBlockConstant(latt_vec2);
|
||
}
|
||
}
|
||
|
||
// Configure solver
|
||
ceres::Solver::Options options;
|
||
options.linear_solver_type = ceres::DENSE_QR;
|
||
options.minimizer_progress_to_stdout = false;
|
||
options.max_solver_time_in_seconds = data.max_time;
|
||
options.logging_type = ceres::LoggingType::SILENT;
|
||
options.num_threads = 1; // Fix threads to 1, as this runs in multi-threaded context
|
||
ceres::Solver::Summary summary;
|
||
|
||
// Run optimization
|
||
ceres::Solve(options, &problem, &summary);
|
||
|
||
if (data.refine_beam_center) {
|
||
data.beam_corr_x = data.geom.GetBeamX_pxl() - beam[0];
|
||
data.beam_corr_y = data.geom.GetBeamY_pxl() - beam[1];
|
||
data.geom.BeamX_pxl(beam[0]).BeamY_pxl(beam[1]);
|
||
}
|
||
|
||
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.axis && data.refine_rotation_axis)
|
||
data.axis.value().Axis(Coord(rot_vec[0], rot_vec[1], rot_vec[2]));
|
||
|
||
if (data.crystal_system == gemmi::CrystalSystem::Orthorhombic)
|
||
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, PI / 2.0, PI / 2.0, PI / 2.0);
|
||
else if (data.crystal_system == gemmi::CrystalSystem::Tetragonal) {
|
||
latt_vec1[1] = latt_vec1[0];
|
||
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, PI / 2.0, PI / 2.0, PI / 2.0);
|
||
} else if (data.crystal_system == gemmi::CrystalSystem::Cubic) {
|
||
latt_vec1[1] = latt_vec1[0];
|
||
latt_vec1[2] = latt_vec1[0];
|
||
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, PI / 2.0, PI / 2.0, PI / 2.0);
|
||
} else if (data.crystal_system == gemmi::CrystalSystem::Hexagonal) {
|
||
latt_vec1[1] = latt_vec1[0];
|
||
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1,PI / 2.0, PI / 2.0, 2.0 * PI / 3.0);
|
||
} else if (data.crystal_system == gemmi::CrystalSystem::Monoclinic) {
|
||
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, PI / 2.0, latt_vec2[0], PI / 2.0);
|
||
} else {
|
||
// Triclinic via the same generic builder
|
||
data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, latt_vec2[0], latt_vec2[1], latt_vec2[2]);
|
||
}
|
||
return true;
|
||
} catch (...) {
|
||
// Convergence problems, likely not updated
|
||
return false;
|
||
}
|
||
}
|
||
|
||
bool XtalOptimizer(XtalOptimizerData &data, const std::vector<std::vector<SpotToSave>> &spots) {
|
||
if (!XtalOptimizerInternal(data, spots, 0.3))
|
||
return false;
|
||
XtalOptimizerInternal(data, spots, 0.2);
|
||
return XtalOptimizerInternal(data, spots, 0.1);
|
||
}
|
||
|
||
bool XtalOptimizerRotationOnly(XtalOptimizerData &data,
|
||
const std::vector<SpotToSave> &spots,
|
||
const float tolerance) {
|
||
try {
|
||
// Parameter: angle-axis for the extra rotation. Identity == {0,0,0}.
|
||
double rot_aa[3] = {0.0, 0.0, 0.0};
|
||
|
||
// Spot selection by current indexing (same approach as XtalOptimizerInternal)
|
||
const Coord a0 = data.latt.Vec0();
|
||
const Coord b0 = data.latt.Vec1();
|
||
const Coord c0 = data.latt.Vec2();
|
||
|
||
const float tol_sq = tolerance * tolerance;
|
||
|
||
ceres::Problem problem;
|
||
|
||
for (const auto &pt : spots) {
|
||
if (!data.index_ice_rings && pt.ice_ring)
|
||
continue;
|
||
|
||
// Compute fractional HKL using the CURRENT lattice
|
||
Coord recip_index = pt.ReciprocalCoord(data.geom);
|
||
if (data.axis.has_value())
|
||
recip_index = data.axis->GetTransformationAngle(pt.phi) * recip_index;
|
||
|
||
const double h_fp = static_cast<double>(recip_index * a0);
|
||
const double k_fp = static_cast<double>(recip_index * b0);
|
||
const double l_fp = static_cast<double>(recip_index * c0);
|
||
|
||
const double h = std::round(h_fp);
|
||
const double k = std::round(k_fp);
|
||
const double l = std::round(l_fp);
|
||
|
||
const double norm_sq =
|
||
(h - h_fp) * (h - h_fp) +
|
||
(k - k_fp) * (k - k_fp) +
|
||
(l - l_fp) * (l - l_fp);
|
||
|
||
if (norm_sq > static_cast<double>(tol_sq))
|
||
continue;
|
||
|
||
// s_obs must be in the same reference frame as the
|
||
// predicted reciprocal vector (h·a* + k·b* + l·c*), which is the
|
||
// phi=0 crystal frame. Apply the same goniometer back-rotation
|
||
// that was used above for the HKL assignment.
|
||
Coord s_obs = data.geom.DetectorToRecip(pt.x, pt.y);
|
||
if (data.axis.has_value())
|
||
s_obs = data.axis->GetTransformationAngle(pt.phi) * s_obs;
|
||
|
||
auto *cost =
|
||
new ceres::AutoDiffCostFunction<XtalResidualRotationOnlyPrecomp, 3, 3>(
|
||
new XtalResidualRotationOnlyPrecomp(s_obs, data.latt, h, k, l)
|
||
);
|
||
|
||
problem.AddResidualBlock(cost, nullptr, rot_aa);
|
||
}
|
||
|
||
if (problem.NumResidualBlocks() < data.min_spots)
|
||
return false;
|
||
|
||
// Regularization: prefer the smallest rotation correction that fits the
|
||
// data. This is essential when spots are nearly coplanar in reciprocal
|
||
// space (e.g. still images), where the rotation component perpendicular
|
||
// to the scattering plane is otherwise underdetermined.
|
||
// The weight is in Å⁻¹ rad⁻¹; tune relative to your typical residual.
|
||
{
|
||
const double reg_weight = 0.05; // e.g. 0.05
|
||
problem.AddResidualBlock(
|
||
new ceres::AutoDiffCostFunction<RotationNormRegularizer, 3, 3>(
|
||
new RotationNormRegularizer(reg_weight)),
|
||
nullptr, rot_aa);
|
||
}
|
||
|
||
ceres::Solver::Options options;
|
||
options.linear_solver_type = ceres::DENSE_QR;
|
||
options.minimizer_progress_to_stdout = false;
|
||
options.max_solver_time_in_seconds = data.max_time;
|
||
options.logging_type = ceres::LoggingType::SILENT;
|
||
options.num_threads = 1;
|
||
|
||
ceres::Solver::Summary summary;
|
||
ceres::Solve(options, &problem, &summary);
|
||
|
||
// Apply rotation to direct-lattice vectors.
|
||
// ceres::AngleAxisToRotationMatrix writes a **row-major** 3×3 matrix,
|
||
// and Eigen's << operator also fills row-by-row, so the assignment
|
||
// below is correct without any transposing.
|
||
//
|
||
// Note: for a pure orthogonal rotation R, R⁻ᵀ = R, so rotating the
|
||
// direct-lattice vectors (A, B, C) by R is exactly equivalent to
|
||
// rotating the reciprocal vectors (a*, b*, c*) by the same R. No
|
||
// transpose or inversion of R is needed here.
|
||
double R_raw[9];
|
||
ceres::AngleAxisToRotationMatrix(rot_aa, R_raw); // row-major 3x3
|
||
|
||
Eigen::Matrix3d R;
|
||
R << R_raw[0], R_raw[3], R_raw[6],
|
||
R_raw[1], R_raw[4], R_raw[7],
|
||
R_raw[2], R_raw[5], R_raw[8];
|
||
|
||
const Eigen::Vector3d A(a0.x, a0.y, a0.z);
|
||
const Eigen::Vector3d B(b0.x, b0.y, b0.z);
|
||
const Eigen::Vector3d C(c0.x, c0.y, c0.z);
|
||
|
||
const Eigen::Vector3d A2 = R * A;
|
||
const Eigen::Vector3d B2 = R * B;
|
||
const Eigen::Vector3d C2 = R * C;
|
||
|
||
data.latt = CrystalLattice(
|
||
Coord(static_cast<float>(A2.x()), static_cast<float>(A2.y()), static_cast<float>(A2.z())),
|
||
Coord(static_cast<float>(B2.x()), static_cast<float>(B2.y()), static_cast<float>(B2.z())),
|
||
Coord(static_cast<float>(C2.x()), static_cast<float>(C2.y()), static_cast<float>(C2.z()))
|
||
);
|
||
|
||
double theta = std::sqrt(rot_aa[0] * rot_aa[0] + rot_aa[1] * rot_aa[1] + rot_aa[2] * rot_aa[2]);
|
||
data.angle_corr = theta;
|
||
if (theta > 1e-6) {
|
||
Coord rot;
|
||
rot.x = rot_aa[0] / theta;
|
||
rot.y = rot_aa[1] / theta;
|
||
rot.z = rot_aa[2] / theta;
|
||
data.angle_axis = rot;
|
||
} else
|
||
data.angle_axis.reset();
|
||
|
||
return true;
|
||
} catch (...) {
|
||
return false;
|
||
}
|
||
} |