diff --git a/image_analysis/geom_refinement/XtalOptimizer.cpp b/image_analysis/geom_refinement/XtalOptimizer.cpp index d4d0c136..1d529aca 100644 --- a/image_analysis/geom_refinement/XtalOptimizer.cpp +++ b/image_analysis/geom_refinement/XtalOptimizer.cpp @@ -138,9 +138,9 @@ struct XtalResidual { const T L1 = e_uc_len[1]; const T L2 = e_uc_len[2]; - T col0_unrot[3] = {L0 * B(0, 0), L1 * B(1, 0), L2 * B(2, 0)}; - T col1_unrot[3] = {L0 * B(0, 1), L1 * B(1, 1), L2 * B(2, 1)}; - T col2_unrot[3] = {L0 * B(0, 2), L1 * B(1, 2), L2 * B(2, 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); @@ -298,7 +298,7 @@ void LatticeToRodriguesAndLengths_Hex(const CrystalLattice &latt, double rod[3], // Extract rotation (Rodrigues), lengths (a,b,c) and beta (rad) for monoclinic (unique axis b). // Frame choice: e2 aligned with b; e1 from a projected orthogonal to e2; e3 = e1 x e2. -inline void LatticeToRodriguesLengthsBeta_Mono(const CrystalLattice &latt, +void LatticeToRodriguesLengthsBeta_Mono(const CrystalLattice &latt, double rod[3], double lengths[3], double &beta_rad) { @@ -325,7 +325,7 @@ inline void LatticeToRodriguesLengthsBeta_Mono(const CrystalLattice &latt, beta_rad = std::acos(cos_beta); // Recover R from the same forward model used in refinement: - // L ≈ R * D(a,b,c) * B(beta) => R ≈ L * (D*B)^-1 + // L ≈ R * B(beta) * D(a,b,c) => R ≈ L * (B*D)^-1 Eigen::Matrix3d L; L.col(0) = A; L.col(1) = Bv; @@ -336,14 +336,13 @@ inline void LatticeToRodriguesLengthsBeta_Mono(const CrystalLattice &latt, Bmono(2, 2) = std::sin(beta_rad); Eigen::DiagonalMatrix D(lengths[0], lengths[1], lengths[2]); - Eigen::Matrix3d M = D * Bmono; + Eigen::Matrix3d M = Bmono * D; Eigen::Matrix3d R_est = Eigen::Matrix3d::Identity(); if (std::abs(M.determinant()) > 1e-15) { R_est = L * M.inverse(); } - // Project to nearest proper rotation (polar decomposition via SVD) Eigen::JacobiSVD svd(R_est, Eigen::ComputeFullU | Eigen::ComputeFullV); Eigen::Matrix3d R = svd.matrixU() * svd.matrixV().transpose(); if (R.determinant() < 0.0) { @@ -373,13 +372,15 @@ CrystalLattice AngleAxisAndLengthsToLattice(const double rod[3], const double le Bhex(1, 1) = sqrt(3) / 2; } - auto latt = R * D * Bhex; + // IMPORTANT: scale columns (a,b,c) by multiplying on the right with D. + auto latt = R * Bhex * 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))); } -inline CrystalLattice AngleAxisLengthsBetaToLattice_Mono(const double rod[3], +CrystalLattice AngleAxisLengthsBetaToLattice_Mono(const double rod[3], const double lengths[3], double beta_rad) { const Eigen::Vector3d r(rod[0], rod[1], rod[2]); @@ -391,11 +392,15 @@ inline CrystalLattice AngleAxisLengthsBetaToLattice_Mono(const double rod[3], const Eigen::DiagonalMatrix D(lengths[0], lengths[1], lengths[2]); 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); - Eigen::Matrix3d latt = R * D * B; + // 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))); @@ -640,7 +645,7 @@ bool XtalOptimizerInternal(XtalOptimizerData &data, B(2, 2) = cz; Eigen::DiagonalMatrix D(latt_vec1[0], latt_vec1[1], latt_vec1[2]); - Eigen::Matrix3d latt = R * D * B; + 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)), @@ -661,3 +666,4 @@ bool XtalOptimizer(XtalOptimizerData &data, const std::vector &spots XtalOptimizerInternal(data, spots, 0.2); return XtalOptimizerInternal(data, spots, 0.1); } + diff --git a/image_analysis/geom_refinement/XtalOptimizer.h b/image_analysis/geom_refinement/XtalOptimizer.h index c78ff96f..5348bca1 100644 --- a/image_analysis/geom_refinement/XtalOptimizer.h +++ b/image_analysis/geom_refinement/XtalOptimizer.h @@ -42,6 +42,14 @@ struct XtalOptimizerData { void LatticeToRodriguesAndLengths_GS(const CrystalLattice &latt, double rod[3], double lengths[3]); void LatticeToRodriguesAndLengths_Hex(const CrystalLattice &latt, double rod[3], double ac[3]); +void LatticeToRodriguesLengthsBeta_Mono(const CrystalLattice &latt, + double rod[3], + 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); diff --git a/tests/XtalOptimizerTest.cpp b/tests/XtalOptimizerTest.cpp index 94e2705c..c14bc846 100644 --- a/tests/XtalOptimizerTest.cpp +++ b/tests/XtalOptimizerTest.cpp @@ -436,7 +436,7 @@ TEST_CASE("XtalOptimizer_monoclinic") { .PoniRot2_rad(0.02) .DetectorDistance_mm(200); - CrystalLattice latt_i(50,60,70,90,96,90); + CrystalLattice latt_i(50,60,70,90,110,90); auto uc_i = latt_i.GetUnitCell(); @@ -455,7 +455,7 @@ TEST_CASE("XtalOptimizer_monoclinic") { } XtalOptimizerData xtal_opt; - xtal_opt.latt = CrystalLattice(49.5, 60.5, 69.8, 90, 95.5, 90); + xtal_opt.latt = CrystalLattice(49.5, 60.5, 69.8, 90, 110, 90); xtal_opt.geom.BeamX_pxl(1007).BeamY_pxl(990).DetectorDistance_mm(200) .PoniRot1_rad(0.01).PoniRot2_rad(0.02); xtal_opt.crystal_system = gemmi::CrystalSystem::Monoclinic; @@ -714,4 +714,153 @@ TEST_CASE("XtalOptimizer_refine_rotation_axis") { CHECK(fabsf(xtal_opt.axis->GetAxis().x - 1.0) < 0.01f); CHECK(fabsf(xtal_opt.axis->GetAxis().y) < 0.01f); CHECK(fabsf(xtal_opt.axis->GetAxis().z) < 0.01f); -} \ No newline at end of file +} + +// --- helpers for lattice sanity tests --- +#include + +namespace { +Eigen::Vector3d to_eigen(const Coord& v) { + return {v[0], v[1], v[2]}; +} + +double angle_rad(const Eigen::Vector3d& a, const Eigen::Vector3d& b) { + const double na = a.norm(); + const double nb = b.norm(); + if (na == 0.0 || nb == 0.0) + return 0.0; + double c = a.dot(b) / (na * nb); + c = std::max(-1.0, std::min(1.0, c)); + return std::acos(c); +} + +// Compare two lattices up to a global rotation: compare Gram matrices G = L^T L (rotation-invariant). +Eigen::Matrix3d gram(const CrystalLattice& latt) { + const Eigen::Vector3d A = to_eigen(latt.Vec0()); + const Eigen::Vector3d B = to_eigen(latt.Vec1()); + const Eigen::Vector3d C = to_eigen(latt.Vec2()); + Eigen::Matrix3d G; + G(0,0) = A.dot(A); G(0,1) = A.dot(B); G(0,2) = A.dot(C); + G(1,0) = B.dot(A); G(1,1) = B.dot(B); G(1,2) = B.dot(C); + G(2,0) = C.dot(A); G(2,1) = C.dot(B); G(2,2) = C.dot(C); + return G; +} + + void check_gram_close(const CrystalLattice& a, + const CrystalLattice& b, + double abs_eps, + double rel_eps) { + const Eigen::Matrix3d Ga = gram(a); + const Eigen::Matrix3d Gb = gram(b); + + for (int r = 0; r < 3; ++r) { + for (int c = 0; c < 3; ++c) { + const double va = Ga(r, c); + const double vb = Gb(r, c); + + // Scale for relative error; avoid blowing up around zero. + const double scale = std::max({1.0, std::abs(va), std::abs(vb)}); + const double tol = std::max(abs_eps, rel_eps * scale); + + INFO("G(" << r << "," << c << ") va=" << va << " vb=" << vb + << " scale=" << scale << " tol=" << tol); + + CHECK(va == Catch::Approx(vb).margin(tol)); + } + } +} +} // namespace + +TEST_CASE("XtalOptimizer Lattice param roundtrip (GS) preserves Gram matrix") { + // Non-orthogonal, irregular basis, but still a valid lattice + CrystalLattice latt_i(Coord(40, 1, 2), + Coord( 3, 50, -4), + Coord(-5, 6, 80)); + + double rod[3]{}, lengths[3]{}; + LatticeToRodriguesAndLengths_GS(latt_i, rod, lengths); + CrystalLattice latt_o = AngleAxisAndLengthsToLattice(rod, lengths, false); + + // This parametrization only keeps "lengths + rotation", i.e. it cannot reproduce shear. + // So we *do not* compare Gram matrices here. + // Instead, sanity-check: reconstructed vectors have the requested lengths. + CHECK(latt_o.Vec0().Length() == Catch::Approx(lengths[0]).margin(1e-9)); + CHECK(latt_o.Vec1().Length() == Catch::Approx(lengths[1]).margin(1e-9)); + CHECK(latt_o.Vec2().Length() == Catch::Approx(lengths[2]).margin(1e-9)); +} + +TEST_CASE("XtalOptimizer Lattice param roundtrip (Hex) preserves unit cell") { + Coord a = Coord(40, 0, 0); + Coord b = Coord(40 / 2.0, 40 * std::sqrt(3) / 2.0, 0); + Coord c = Coord(0, 0, 70); + + // Apply an arbitrary rotation to ensure the rod extraction is meaningful + RotMatrix R(1.0, Coord(0, 1, 1)); + CrystalLattice latt_i(R * a, R * b, R * c); + + double rod[3]{}, ac[3]{}; + LatticeToRodriguesAndLengths_Hex(latt_i, rod, ac); + CrystalLattice latt_o = AngleAxisAndLengthsToLattice(rod, ac, true); + + auto uc_o = latt_o.GetUnitCell(); + CHECK(uc_o.a == Catch::Approx(40.0).margin(1e-6)); + CHECK(uc_o.b == Catch::Approx(40.0).margin(1e-6)); + CHECK(uc_o.c == Catch::Approx(70.0).margin(1e-6)); + CHECK(uc_o.alpha == Catch::Approx(90.0).margin(1e-6)); + CHECK(uc_o.beta == Catch::Approx(90.0).margin(1e-6)); + CHECK(uc_o.gamma == Catch::Approx(120.0).margin(1e-6)); +} + +TEST_CASE("XtalOptimizer Monoclinic param roundtrip preserves Gram matrix (beta far from 90)") { + struct Case { double beta_deg; }; + const std::vector cases = { + {60.0}, {75.0}, {115.0}, {130.0} + }; + + for (const auto& cs : cases) { + INFO("beta_deg=" << cs.beta_deg); + + // Start from a clean monoclinic cell in its conventional setting (unique axis b). + CrystalLattice latt0(50, 60, 70, 90, cs.beta_deg, 90); + + // Now apply a TRUE global rotation: rotate each basis vector (left-multiply). + RotMatrix R(0.7, Coord(0.3, 0.9, 0.1)); + CrystalLattice latt_i(R * latt0.Vec0(), + R * latt0.Vec1(), + R * latt0.Vec2()); + + double rod[3]{}, lengths[3]{}, beta_rad = 0.0; + LatticeToRodriguesLengthsBeta_Mono(latt_i, rod, lengths, beta_rad); + + // Basic sanity + CHECK(lengths[0] == Catch::Approx(50.0).margin(1e-6)); + CHECK(lengths[1] == Catch::Approx(60.0).margin(1e-6)); + 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); + + // Rotation-invariant check: Gram matrices match. + check_gram_close(latt_i, latt_o, /*abs_eps=*/5e-4, /*rel_eps=*/1e-10); + + // Also check the unit-cell angles we expect for monoclinic(unique b). + auto uc_o = latt_o.GetUnitCell(); + CHECK(uc_o.alpha == Catch::Approx(90.0).margin(1e-4)); + CHECK(uc_o.gamma == Catch::Approx(90.0).margin(1e-4)); + CHECK(uc_o.beta == Catch::Approx(cs.beta_deg).margin(1e-4)); + } +} + +TEST_CASE("XtalOptimizer Monoclinic beta geometry: extracted beta equals angle(a,c)") { + // This isolates only the beta definition, independent of other choices. + CrystalLattice latt_i(50, 60, 70, 90, 130, 90); + + const Eigen::Vector3d A = to_eigen(latt_i.Vec0()); + const Eigen::Vector3d C = to_eigen(latt_i.Vec2()); + const double beta_geom = angle_rad(A, C); + + double rod[3]{}, lengths[3]{}, beta_rad = 0.0; + LatticeToRodriguesLengthsBeta_Mono(latt_i, rod, lengths, beta_rad); + + CHECK(beta_rad == Catch::Approx(beta_geom).margin(1e-12)); +}