From ecdb7048a018a2dbf33bf27088045401ddcbb728 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Mon, 15 Jun 2026 14:11:26 +0200 Subject: [PATCH] PixelRefine: drop crystal-system idealisation, use the indexed cell as-is 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 --- image_analysis/IndexAndRefine.cpp | 3 - .../pixel_refinement/PixelRefine.cpp | 153 +++++------------- image_analysis/pixel_refinement/PixelRefine.h | 7 +- 3 files changed, 41 insertions(+), 122 deletions(-) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 4800e8bb..ebe85b0c 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -439,9 +439,6 @@ bool IndexAndRefine::PixelRefineIntegrate(DataMessage &msg, PixelRefineData prd; prd.geom = outcome.experiment.GetDiffractionGeometry(); prd.latt = *outcome.lattice_candidate; - prd.crystal_system = outcome.symmetry.crystal_system; - if (prd.crystal_system == gemmi::CrystalSystem::Trigonal) - prd.crystal_system = gemmi::CrystalSystem::Hexagonal; prd.centering = outcome.symmetry.centering; if (const auto bw = experiment.GetBandwidthFWHM()) prd.bandwidth = bw.value() / 2.3548; // FWHM -> sigma diff --git a/image_analysis/pixel_refinement/PixelRefine.cpp b/image_analysis/pixel_refinement/PixelRefine.cpp index 4f45c7aa..b9014d82 100644 --- a/image_analysis/pixel_refinement/PixelRefine.cpp +++ b/image_analysis/pixel_refinement/PixelRefine.cpp @@ -11,9 +11,6 @@ #include #include -#include - -#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 -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 &e_pred_recip, - Eigen::Matrix &n_radial, T &q_sq) { - Eigen::Matrix e_uc_len = Eigen::Matrix::Zero(); - Eigen::Matrix Bmat = Eigen::Matrix::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 A(col0_rot[0], col0_rot[1], col0_rot[2]); - const Eigen::Matrix Bv(col1_rot[0], col1_rot[1], col1_rot[2]); - const Eigen::Matrix C(col2_rot[0], col2_rot[1], col2_rot[2]); - - const Eigen::Matrix BxC = Bv.cross(C); - const Eigen::Matrix CxA = C.cross(A); - const Eigen::Matrix 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 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(beam, &dist_mm, detector_rot, obs_x, obs_y, pixel_size, inv_lambda, e_obs); Eigen::Vector3d e_pred, n_radial; - if (!PredictedNode(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(x), static_cast(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; diff --git a/image_analysis/pixel_refinement/PixelRefine.h b/image_analysis/pixel_refinement/PixelRefine.h index 366d8ca3..b180e1de 100644 --- a/image_analysis/pixel_refinement/PixelRefine.h +++ b/image_analysis/pixel_refinement/PixelRefine.h @@ -58,7 +58,6 @@ struct PixelRefineData { // --- model state (input as initial guess, output as refined result) --- DiffractionGeometry geom; // fixed (refined upstream by XtalOptimizer) CrystalLattice latt; // fixed - gemmi::CrystalSystem crystal_system = gemmi::CrystalSystem::Triclinic; char centering = 'P'; double B_factor = 0.0; // Debye-Waller B (A^2), refined @@ -113,9 +112,9 @@ class PixelRefine { const HKLKeyGenerator hkl_key_generator; std::map reference_data; - // Fills the fixed geometry + symmetry-aware lattice parametrization (beam, - // distance, detector tilt, and the Rodrigues orientation / cell-length / angle - // vectors) from the current model state, for the per-pixel geometry evaluation. + // Fills the fixed geometry (beam, distance, detector tilt) and the three + // real-space lattice column vectors from the current model state, for the + // per-pixel geometry evaluation. void BuildParameterBlocks(const PixelRefineData &data, double beam[2], double &dist_mm, double detector_rot[2],