XtalOptimizer: Have one function to rebuild lattice from three arrays

image_analysis/geom_refinement/XtalOptimizer.cpp

@@ -1,6 +1,8 @@
 // SPDX-FileCopyrightText: 2025 Filip Leonarski, Paul Scherrer Institute <filip.leonarski@psi.ch>
 // SPDX-License-Identifier: GPL-3.0-only
 
+#include <Eigen/Dense>
+
 #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<double, 3> 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<double, 3> 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<double, 3> 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;

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<SpotToSave> &spots);
 

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);