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

PixelRefine: drop crystal-system idealisation, use the indexed cell as-is
Build Packages / build:rpm (ubuntu2404_nocuda) (push) Successful in 22m4s
Details
Build Packages / build:rpm (rocky8_nocuda) (push) Successful in 24m6s
Details
Build Packages / build:rpm (rocky9_nocuda) (push) Successful in 26m1s
Details
Build Packages / build:rpm (ubuntu2204_nocuda) (push) Successful in 25m21s
Details
Build Packages / build:rpm (rocky8_sls9) (push) Successful in 25m21s
Details
Build Packages / build:rpm (rocky9_sls9) (push) Successful in 25m16s
Details
Build Packages / build:rpm (rocky8) (push) Successful in 24m23s
Details
Build Packages / build:rpm (ubuntu2404) (push) Successful in 22m48s
Details
Build Packages / build:rpm (ubuntu2204) (push) Successful in 25m18s
Details
Build Packages / build:rpm (rocky9) (push) Successful in 27m54s
Details
Build Packages / Generate python client (push) Successful in 32s
Details
Build Packages / Create release (push) Skipped
Details
Build Packages / Build documentation (push) Successful in 1m22s
Details
Build Packages / XDS test (durin plugin) (push) Successful in 19m3s
Details
Build Packages / DIALS test (push) Successful in 30m12s
Details
Build Packages / XDS test (neggia plugin) (push) Successful in 13m46s
Details
Build Packages / XDS test (JFJoch plugin) (push) Successful in 15m34s
Details
Build Packages / Unit tests (push) Successful in 2h33m4s
Details

Browse Source 
PixelRefine no longer refines the cell, so the symmetry-constrained cell
parameterisation inherited from XtalOptimizer has no manifold to constrain.
The decompose->rebuild round-trip (Gram-Schmidt orientation + symmetry B
matrix) merely reconstructed the lattice columns it started from, and for
non-triclinic systems it re-idealised the indexed cell (averaging a,b for
tetragonal; a,b,c for cubic; forcing alpha=gamma=90).

Replace both six-way switches (PredictedNode and BuildParameterBlocks) with a
single path: take the three real-space lattice columns (latt.Vec0/1/2())
directly and form the reciprocal node via the general cross-product inverse.
This reproduces every crystal system exactly from the actual cell, is more
faithful (no re-idealisation), and removes the crystal_system field plus two
now-unused includes. PredictedNode de-templated (only ever called with double).
Crystal symmetry still lives where it belongs: indexing (upstream) and merging
(downstream via the space-group HKL key).

A/B (both lyso P4_3 2_1 2 crystals, refined under tetragonal constraint, so the
old idealisation was already a no-op): stat tables bit-identical across all
shells -- crystal 2 CC1/2 94.6% / CCref 92.5%, jet CC1/2 91.9% / CCref 55.8% --
the only delta is 4th-digit float-ordering noise in the error-model b
coefficient (0.8143->0.8141, same ISa). Same merged intensities => R-free
unchanged. Net -81 lines.

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-06-15 14:11:26 +02:00
co-authored by Claude Opus 4.8
parent 8a9d80eb71
commit ecdb7048a0
3 changed files with 41 additions and 122 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/pixel_refinement/PixelRefine.cpp
+38 -115
View File
@@ -11,9 +11,6 @@
#include <Eigen/Dense>
#include <ceres/ceres.h>
#include <ceres/rotation.h>
#include "../geom_refinement/LatticeReduction.h"
namespace {
@@ -186,89 +183,36 @@ void ObservedRecip(const T *beam, const T *distance_mm, const T *detector_rot,
}
// Per-reflection: predicted node g_hkl, |g_hkl|^2, and the Ewald-sphere normal.
// This is the expensive part (symmetry-aware B matrix, three rotations, cross
// products) - it depends only on the lattice (p0,p1,p2) and hkl, so for a whole
// shoebox it can be computed once. Convention identical to XtalOptimizer.
template<typename T>
bool PredictedNode(const T *p0, const T *p1, const T *p2,
double exp_h, double exp_k, double exp_l,
gemmi::CrystalSystem symmetry, double inv_lambda,
Eigen::Matrix<T, 3, 1> &e_pred_recip,
Eigen::Matrix<T, 3, 1> &n_radial, T &q_sq) {
Eigen::Matrix<T, 3, 1> e_uc_len = Eigen::Matrix<T, 3, 1>::Zero();
Eigen::Matrix<T, 3, 3> Bmat = Eigen::Matrix<T, 3, 3>::Identity();
// (a0,a1,a2) are the three real-space lattice column vectors in the lab frame
// (latt.Vec0/1/2()), used exactly as indexed: PixelRefine does not refine the cell,
// so there is no symmetry manifold to constrain - a general (triclinic) reciprocal
// inverse reproduces every crystal system from the actual cell. Depends only on the
// lattice and hkl, so for a whole shoebox it is computed once.
bool PredictedNode(const double *a0, const double *a1, const double *a2,
double exp_h, double exp_k, double exp_l, double inv_lambda,
Eigen::Vector3d &e_pred_recip,
Eigen::Vector3d &n_radial, double &q_sq) {
const Eigen::Vector3d A(a0[0], a0[1], a0[2]);
const Eigen::Vector3d Bv(a1[0], a1[1], a1[2]);
const Eigen::Vector3d C(a2[0], a2[1], a2[2]);
if (symmetry == gemmi::CrystalSystem::Hexagonal) {
e_uc_len << p1[0], p1[0], p1[2];
Bmat(0, 1) = T(-0.5);
Bmat(1, 1) = T(sqrt(3.0) / 2.0);
} 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) {
e_uc_len << p1[0], p1[1], p1[2];
Bmat(0, 2) = ceres::cos(p2[0]);
Bmat(2, 2) = ceres::sin(p2[0]);
} else {
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]);
const Eigen::Vector3d BxC = Bv.cross(C);
const Eigen::Vector3d CxA = C.cross(A);
const Eigen::Vector3d AxB = A.cross(Bv);
e_uc_len << p1[0], p1[1], p1[2];
Bmat(0, 1) = cg;
Bmat(1, 1) = sg;
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);
Bmat(0, 2) = cx;
Bmat(1, 2) = cy;
Bmat(2, 2) = cz;
}
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] = {Bmat(0, 0) * L0, Bmat(1, 0) * L0, Bmat(2, 0) * L0};
T col1_unrot[3] = {Bmat(0, 1) * L1, Bmat(1, 1) * L1, Bmat(2, 1) * L1};
T col2_unrot[3] = {Bmat(0, 2) * L2, Bmat(1, 2) * L2, Bmat(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 Vol = A.dot(BxC);
if (ceres::abs(Vol) < T(1e-12))
const double Vol = A.dot(BxC);
if (std::abs(Vol) < 1e-12)
return false;
const T invV = T(1) / Vol;
const double invV = 1.0 / Vol;
e_pred_recip = (BxC * T(exp_h) + CxA * T(exp_k) + AxB * T(exp_l)) * invV;
e_pred_recip = (BxC * exp_h + CxA * exp_k + AxB * exp_l) * invV;
q_sq = e_pred_recip.squaredNorm();
// Ewald sphere centre at -k_i = (0,0,-inv_lambda); radial normal at g_hkl.
const Eigen::Matrix<T, 3, 1> S_pred(
e_pred_recip[0],
e_pred_recip[1],
e_pred_recip[2] + T(inv_lambda));
const T S_pred_norm = S_pred.norm();
if (S_pred_norm < T(1e-10))
const Eigen::Vector3d S_pred(
e_pred_recip[0], e_pred_recip[1], e_pred_recip[2] + inv_lambda);
const double S_pred_norm = S_pred.norm();
if (S_pred_norm < 1e-10)
return false;
n_radial = S_pred / S_pred_norm;
@@ -281,7 +225,7 @@ bool PredictedNode(const T *p0, const T *p1, const T *p2,
// eps_radial = deviation along the Ewald normal (the partiality direction)
// eps_tang_sq = squared deviation in the tangent plane (the profile direction)
bool GeometryProbe(double obs_x, double obs_y, double lambda, double pixel_size,
int h, int k, int l, gemmi::CrystalSystem symmetry,
int h, int k, int l,
const double beam[2], double dist_mm, const double detector_rot[2],
const double p0[3], const double p1[3], const double p2[3],
double &q_sq, double &eps_radial, double &eps_tang_sq) {
@@ -290,7 +234,7 @@ bool GeometryProbe(double obs_x, double obs_y, double lambda, double pixel_size,
ObservedRecip<double>(beam, &dist_mm, detector_rot, obs_x, obs_y, pixel_size, inv_lambda, e_obs);
Eigen::Vector3d e_pred, n_radial;
if (!PredictedNode<double>(p0, p1, p2, h, k, l, symmetry, inv_lambda, e_pred, n_radial, q_sq))
if (!PredictedNode(p0, p1, p2, h, k, l, inv_lambda, e_pred, n_radial, q_sq))
return false;
const Eigen::Vector3d delta_q = e_obs - e_pred;
@@ -359,37 +303,16 @@ void PixelRefine::BuildParameterBlocks(const PixelRefineData &data,
detector_rot[0] = data.geom.GetPoniRot1_rad();
detector_rot[1] = data.geom.GetPoniRot2_rad();
for (int i = 0; i < 3; ++i)
latt_vec0[i] = latt_vec1[i] = latt_vec2[i] = 0.0;
double beta = data.latt.GetUnitCell().beta;
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;
break;
default: {
LatticeToRodriguesAndLengths_GS(data.latt, latt_vec0, latt_vec1);
const auto uc = data.latt.GetUnitCell();
latt_vec2[0] = uc.alpha * M_PI / 180.0;
latt_vec2[1] = uc.beta * M_PI / 180.0;
latt_vec2[2] = uc.gamma * M_PI / 180.0;
break;
}
// The three real-space lattice columns in the lab frame, used as-is. PixelRefine
// does not refine the cell, so there is no symmetry manifold to constrain; the
// indexed cell is taken exactly, with no re-idealisation to ideal symmetry.
const Coord &a = data.latt.Vec0();
const Coord &b = data.latt.Vec1();
const Coord &c = data.latt.Vec2();
for (int i = 0; i < 3; ++i) {
latt_vec0[i] = a[i];
latt_vec1[i] = b[i];
latt_vec2[i] = c[i];
}
}
@@ -532,7 +455,7 @@ void PixelRefine::Run(const T *image,
double sw = 0.0, sw_et2 = 0.0;
for (const auto &px : g.pixels) {
double q_sq, eps_r, eps_t_sq;
if (!GeometryProbe(px.x, px.y, lambda, pixel_size, g.h, g.k, g.l, data.crystal_system,
if (!GeometryProbe(px.x, px.y, lambda, pixel_size, g.h, g.k, g.l,
beam, dist_mm, detector_rot, latt_vec0, latt_vec1, latt_vec2,
q_sq, eps_r, eps_t_sq))
continue;
@@ -572,7 +495,7 @@ void PixelRefine::Run(const T *image,
pt_sig.reserve(g.pixels.size());
for (const auto &px : g.pixels) {
double q_sq, eps_r, eps_t_sq;
if (!GeometryProbe(px.x, px.y, lambda, pixel_size, g.h, g.k, g.l, data.crystal_system,
if (!GeometryProbe(px.x, px.y, lambda, pixel_size, g.h, g.k, g.l,
beam, dist_mm, detector_rot, latt_vec0, latt_vec1, latt_vec2,
q_sq, eps_r, eps_t_sq))
continue;
@@ -664,7 +587,7 @@ void PixelRefine::Run(const T *image,
for (int x = cx - radius; x <= cx + radius; ++x) {
double q_sq, eps_r, eps_t_sq;
if (!GeometryProbe(static_cast<double>(x), static_cast<double>(y),
lambda, pixel_size, g.h, g.k, g.l, data.crystal_system,
lambda, pixel_size, g.h, g.k, g.l,
beam, dist_mm, detector_rot, latt_vec0, latt_vec1, latt_vec2,
q_sq, eps_r, eps_t_sq))
continue;