diff --git a/image_analysis/geom_refinement/XtalOptimizer.cpp b/image_analysis/geom_refinement/XtalOptimizer.cpp index 1d529aca..3809d314 100644 --- a/image_analysis/geom_refinement/XtalOptimizer.cpp +++ b/image_analysis/geom_refinement/XtalOptimizer.cpp @@ -1,6 +1,8 @@ // SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute // SPDX-License-Identifier: GPL-3.0-only +#include + #include "XtalOptimizer.h" #include "ceres/ceres.h" #include "ceres/rotation.h" @@ -132,7 +134,7 @@ struct XtalResidual { B(2, 2) = cz; } - // Build unrotated direct lattice columns: (D * B), then rotate them by p0. + // 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]; @@ -358,52 +360,67 @@ void LatticeToRodriguesLengthsBeta_Mono(const CrystalLattice &latt, rod[2] = r.z(); } -CrystalLattice AngleAxisAndLengthsToLattice(const double rod[3], const double lengths[3], bool hex) { +static inline Eigen::Matrix3d B_from_angles(double alpha_rad, double beta_rad, double gamma_rad) { + const double ca = std::cos(alpha_rad); + const double cb = std::cos(beta_rad); + const double cg = std::cos(gamma_rad); + const double sg = std::sin(gamma_rad); + + Eigen::Matrix3d B = Eigen::Matrix3d::Identity(); + + // 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; + + // c vector components (standard crystallography construction) + 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)); + + B(0, 2) = cx; + B(1, 2) = cy; + B(2, 2) = cz; + + return B; +} + +CrystalLattice AngleAxisAndCellToLattice(const double rod[3], + const double lengths[3], + double alpha_rad, + double beta_rad, + double gamma_rad) { const Eigen::Vector3d r(rod[0], rod[1], rod[2]); const double angle = r.norm(); + Eigen::Matrix3d R = Eigen::Matrix3d::Identity(); if (angle > 0.0) R = Eigen::AngleAxisd(angle, r / angle).toRotationMatrix(); + const Eigen::DiagonalMatrix D(lengths[0], lengths[1], lengths[2]); + const Eigen::Matrix3d B = B_from_angles(alpha_rad, beta_rad, gamma_rad); - Eigen::Matrix3d Bhex = Eigen::Matrix3d::Identity(); - if (hex) { - Bhex(0, 1) = -1 / 2.0; - Bhex(1, 1) = sqrt(3) / 2; - } - - // IMPORTANT: scale columns (a,b,c) by multiplying on the right with D. - auto latt = R * Bhex * D; + // IMPORTANT convention: L = R * B * D (scale columns by lengths) + const Eigen::Matrix3d latt = R * B * D; return CrystalLattice(Coord(latt(0, 0), latt(1, 0), latt(2, 0)), Coord(latt(0, 1), latt(1, 1), latt(2, 1)), Coord(latt(0, 2), latt(1, 2), latt(2, 2))); } -CrystalLattice AngleAxisLengthsBetaToLattice_Mono(const double rod[3], - const double lengths[3], - double beta_rad) { - const Eigen::Vector3d r(rod[0], rod[1], rod[2]); - const double angle = r.norm(); - Eigen::Matrix3d R = Eigen::Matrix3d::Identity(); - if (angle > 0.0) - R = Eigen::AngleAxisd(angle, r / angle).toRotationMatrix(); - const Eigen::DiagonalMatrix D(lengths[0], lengths[1], lengths[2]); +CrystalLattice AngleAxisAndLengthsToLattice(const double rod[3], const double lengths[3], bool hex) { + if (!hex) { + return AngleAxisAndCellToLattice(rod, lengths, + /*alpha=*/M_PI / 2.0, + /*beta =*/M_PI / 2.0, + /*gamma=*/M_PI / 2.0); + } - Eigen::Matrix3d B = Eigen::Matrix3d::Identity(); - // Columns are unit directions of (a,b,c) in the unrotated crystal frame. - // For monoclinic (unique b): a along x, b along y, c in x-z with angle beta to a. - // Bmono = [[1,0,cosβ],[0,1,0],[0,0,sinβ]] - B(0, 2) = std::cos(beta_rad); - B(2, 2) = std::sin(beta_rad); - - // IMPORTANT: scale columns (a,b,c) by multiplying on the right with D. - Eigen::Matrix3d latt = R * B * D; - - return CrystalLattice(Coord(latt(0, 0), latt(1, 0), latt(2, 0)), - Coord(latt(0, 1), latt(1, 1), latt(2, 1)), - Coord(latt(0, 2), latt(1, 2), latt(2, 2))); + // Hexagonal: caller must already enforce a=b in `lengths`. + return AngleAxisAndCellToLattice(rod, lengths, + /*alpha=*/M_PI / 2.0, + /*beta =*/M_PI / 2.0, + /*gamma=*/2.0 * M_PI / 3.0); } bool XtalOptimizerInternal(XtalOptimizerData &data, @@ -601,56 +618,24 @@ bool XtalOptimizerInternal(XtalOptimizerData &data, data.axis.value().Axis(Coord(rot_vec[0], rot_vec[1], rot_vec[2])); if (data.crystal_system == gemmi::CrystalSystem::Orthorhombic) - data.latt = AngleAxisAndLengthsToLattice(latt_vec0, latt_vec1, false); + data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, M_PI / 2.0, M_PI / 2.0, M_PI / 2.0); else if (data.crystal_system == gemmi::CrystalSystem::Tetragonal) { latt_vec1[1] = latt_vec1[0]; - data.latt = AngleAxisAndLengthsToLattice(latt_vec0, latt_vec1, false); + data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, M_PI / 2.0, M_PI / 2.0, M_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 = AngleAxisAndLengthsToLattice(latt_vec0, latt_vec1, false); + data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, M_PI / 2.0, M_PI / 2.0, M_PI / 2.0); } else if (data.crystal_system == gemmi::CrystalSystem::Hexagonal) { latt_vec1[1] = latt_vec1[0]; - data.latt = AngleAxisAndLengthsToLattice(latt_vec0, latt_vec1, true); + data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1,M_PI / 2.0, M_PI / 2.0, 2.0 * M_PI / 3.0); } else if (data.crystal_system == gemmi::CrystalSystem::Monoclinic) { - data.latt = AngleAxisLengthsBetaToLattice_Mono(latt_vec0, latt_vec1, latt_vec2[0]); + data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, M_PI / 2.0, latt_vec2[0], M_PI / 2.0); } else { - // Triclinic: reconstruct with generic B from α,β,γ - const Eigen::Vector3d r(latt_vec0[0], latt_vec0[1], latt_vec0[2]); - const double angle = r.norm(); - Eigen::Matrix3d R = Eigen::Matrix3d::Identity(); - if (angle > 0.0) - R = Eigen::AngleAxisd(angle, r / angle).toRotationMatrix(); - - Eigen::Matrix3d B = Eigen::Matrix3d::Identity(); - const double ca = std::cos(latt_vec2[0]); - const double cb = std::cos(latt_vec2[1]); - const double cg = std::cos(latt_vec2[2]); - 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; - - 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)); - - B(0, 2) = cx; - B(1, 2) = cy; - B(2, 2) = cz; - - Eigen::DiagonalMatrix D(latt_vec1[0], latt_vec1[1], latt_vec1[2]); - Eigen::Matrix3d latt = R * B * D; - - data.latt = CrystalLattice(Coord(latt(0, 0), latt(1, 0), latt(2, 0)), - Coord(latt(0, 1), latt(1, 1), latt(2, 1)), - Coord(latt(0, 2), latt(1, 2), latt(2, 2))); + // Triclinic via the same generic builder + data.latt = AngleAxisAndCellToLattice(latt_vec0, latt_vec1, latt_vec2[0], latt_vec2[1], latt_vec2[2]); } + data.latt.Regularize(data.crystal_system); return true; diff --git a/image_analysis/geom_refinement/XtalOptimizer.h b/image_analysis/geom_refinement/XtalOptimizer.h index 5348bca1..79422e6e 100644 --- a/image_analysis/geom_refinement/XtalOptimizer.h +++ b/image_analysis/geom_refinement/XtalOptimizer.h @@ -47,11 +47,11 @@ void LatticeToRodriguesLengthsBeta_Mono(const CrystalLattice &latt, double lengths[3], double &beta_rad); -CrystalLattice AngleAxisLengthsBetaToLattice_Mono(const double rod[3], - const double lengths[3], - double beta_rad); - -CrystalLattice AngleAxisAndLengthsToLattice(const double rod[3], const double lengths[3], bool hex = false); +CrystalLattice AngleAxisAndCellToLattice(const double rod[3], + const double lengths[3], + double alpha_rad, + double beta_rad, + double gamma_rad); bool XtalOptimizer(XtalOptimizerData &data, const std::vector &spots); diff --git a/tests/XtalOptimizerTest.cpp b/tests/XtalOptimizerTest.cpp index 35b951c3..2403ff48 100644 --- a/tests/XtalOptimizerTest.cpp +++ b/tests/XtalOptimizerTest.cpp @@ -494,7 +494,7 @@ TEST_CASE("LatticeToRodrigues") { CHECK(fabs(rod[1]) < 0.001); CHECK(fabs(rod[2]) < 0.001); - auto latt_o = AngleAxisAndLengthsToLattice(rod, lengths, false); + auto latt_o = AngleAxisAndCellToLattice(rod, lengths, M_PI/2, M_PI/2, M_PI/2); CHECK(latt_o.Vec0().Length() == Catch::Approx(40.0)); CHECK(latt_o.Vec1().Length() == Catch::Approx(50.0)); @@ -514,8 +514,7 @@ TEST_CASE("LatticeToRodrigues_irregular") { CHECK(lengths[1] == Catch::Approx(50.0)); CHECK(lengths[2] == Catch::Approx(80.0)); - - auto latt_o = AngleAxisAndLengthsToLattice(rod, lengths, false); + auto latt_o = AngleAxisAndCellToLattice(rod, lengths, M_PI/2, M_PI/2, M_PI/2); CHECK(latt_o.Vec0().Length() == Catch::Approx(40.0)); CHECK(latt_o.Vec1().Length() == Catch::Approx(50.0)); @@ -538,7 +537,7 @@ TEST_CASE("LatticeToRodrigues_Hex") { CHECK(lengths[1] == Catch::Approx(40.0)); CHECK(lengths[2] == Catch::Approx(70.0)); - auto latt_o = AngleAxisAndLengthsToLattice(rod, lengths, true); + auto latt_o = AngleAxisAndCellToLattice(rod, lengths,M_PI / 2.0, M_PI / 2.0, 2.0 * M_PI / 3.0); auto uc_o = latt_o.GetUnitCell(); CHECK(uc_o.a == Catch::Approx(40.0)); CHECK(uc_o.b == Catch::Approx(40.0)); @@ -782,7 +781,7 @@ TEST_CASE("XtalOptimizer Lattice param roundtrip (GS) preserves Gram matrix") { double rod[3]{}, lengths[3]{}; LatticeToRodriguesAndLengths_GS(latt_i, rod, lengths); - CrystalLattice latt_o = AngleAxisAndLengthsToLattice(rod, lengths, false); + auto latt_o = AngleAxisAndCellToLattice(rod, lengths, M_PI/2, M_PI/2, M_PI/2); // This parametrization only keeps "lengths + rotation", i.e. it cannot reproduce shear. // So we *do not* compare Gram matrices here. @@ -803,7 +802,7 @@ TEST_CASE("XtalOptimizer Lattice param roundtrip (Hex) preserves unit cell") { double rod[3]{}, ac[3]{}; LatticeToRodriguesAndLengths_Hex(latt_i, rod, ac); - CrystalLattice latt_o = AngleAxisAndLengthsToLattice(rod, ac, true); + auto latt_o = AngleAxisAndCellToLattice(rod, ac, M_PI/2, M_PI/2, 2*M_PI/3); auto uc_o = latt_o.GetUnitCell(); CHECK(uc_o.a == Catch::Approx(40.0).margin(1e-6)); @@ -841,7 +840,7 @@ TEST_CASE("XtalOptimizer Monoclinic param roundtrip preserves Gram matrix (beta CHECK(lengths[2] == Catch::Approx(70.0).margin(1e-6)); CHECK(beta_rad * 180.0 / M_PI == Catch::Approx(cs.beta_deg).margin(1e-6)); - CrystalLattice latt_o = AngleAxisLengthsBetaToLattice_Mono(rod, lengths, beta_rad); + auto latt_o = AngleAxisAndCellToLattice(rod, lengths, M_PI/2, beta_rad, M_PI/2); // Rotation-invariant check: Gram matrices match. check_gram_close(latt_i, latt_o, /*abs_eps=*/5e-4, /*rel_eps=*/1e-10);