From 7634cd819a742306252bc8d3b2a22c1c692a76fb Mon Sep 17 00:00:00 2001
From: Filip Leonarski <filip.leonarski@psi.ch>
Date: Wed, 18 Feb 2026 20:44:58 +0100
Subject: [PATCH] XtalOptimizer: Fix how monoclinic is generated

---
 .../geom_refinement/XtalOptimizer.cpp         | 53 +++++++++----------
 1 file changed, 26 insertions(+), 27 deletions(-)

diff --git a/image_analysis/geom_refinement/XtalOptimizer.cpp b/image_analysis/geom_refinement/XtalOptimizer.cpp
index 4b05bca4..d4d0c136 100644
--- a/image_analysis/geom_refinement/XtalOptimizer.cpp
+++ b/image_analysis/geom_refinement/XtalOptimizer.cpp
@@ -307,11 +307,11 @@ inline void LatticeToRodriguesLengthsBeta_Mono(const CrystalLattice &latt,
     const Coord c = latt.Vec2();
 
     Eigen::Vector3d A(a[0], a[1], a[2]);
-    Eigen::Vector3d B(b[0], b[1], b[2]);
+    Eigen::Vector3d Bv(b[0], b[1], b[2]);
     Eigen::Vector3d C(c[0], c[1], c[2]);
 
     const double a_len = A.norm();
-    const double b_len = B.norm();
+    const double b_len = Bv.norm();
     const double c_len = C.norm();
 
     lengths[0] = a_len;
@@ -323,35 +323,34 @@ inline void LatticeToRodriguesLengthsBeta_Mono(const CrystalLattice &latt,
     if (a_len > 0.0 && c_len > 0.0)
         cos_beta = std::max(-1.0, std::min(1.0, A.dot(C) / (a_len * c_len)));
     beta_rad = std::acos(cos_beta);
-    Eigen::Vector3d e2, ax;
 
-    // e2 along b (unique axis)
-    if (b_len > 0.0)
-        e2 = B / b_len;
-    else
-        e2 = Eigen::Vector3d::UnitY();
+    // Recover R from the same forward model used in refinement:
+    // L ≈ R * D(a,b,c) * B(beta)  =>  R ≈ L * (D*B)^-1
+    Eigen::Matrix3d L;
+    L.col(0) = A;
+    L.col(1) = Bv;
+    L.col(2) = C;
 
-    if (a_len > 0.0)
-        ax = A / a_len;
-    else
-        ax = Eigen::Vector3d::UnitX();
+    Eigen::Matrix3d Bmono = Eigen::Matrix3d::Identity();
+    Bmono(0, 2) = std::cos(beta_rad);
+    Bmono(2, 2) = std::sin(beta_rad);
 
-    Eigen::Vector3d e1 = ax - (ax.dot(e2)) * e2;
-    double n1 = e1.norm();
-    if (n1 < 1e-15) {
-        // Fallback: use any perpendicular to e2
-        e1 = (std::abs(e2.x()) < 0.9
-                  ? Eigen::Vector3d::UnitX().cross(e2)
-                  : Eigen::Vector3d::UnitY().cross(e2));
+    Eigen::DiagonalMatrix<double, 3> D(lengths[0], lengths[1], lengths[2]);
+    Eigen::Matrix3d M = D * Bmono;
+
+    Eigen::Matrix3d R_est = Eigen::Matrix3d::Identity();
+    if (std::abs(M.determinant()) > 1e-15) {
+        R_est = L * M.inverse();
     }
-    e1.normalize();
-    // e3 completes the right-handed frame
-    Eigen::Vector3d e3 = e1.cross(e2).normalized();
 
-    Eigen::Matrix3d R;
-    R.col(0) = e1;
-    R.col(1) = e2;
-    R.col(2) = e3;
+    // Project to nearest proper rotation (polar decomposition via SVD)
+    Eigen::JacobiSVD<Eigen::Matrix3d> svd(R_est, Eigen::ComputeFullU | Eigen::ComputeFullV);
+    Eigen::Matrix3d R = svd.matrixU() * svd.matrixV().transpose();
+    if (R.determinant() < 0.0) {
+        Eigen::Matrix3d U = svd.matrixU();
+        U.col(2) *= -1.0;
+        R = U * svd.matrixV().transpose();
+    }
 
     Eigen::AngleAxisd aa(R);
     Eigen::Vector3d r = aa.angle() * aa.axis();
@@ -520,7 +519,7 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
         else {
             const double dist_range = 0.1;
             problem.SetParameterLowerBound(&distance_mm, 0, distance_mm * (1.0 - dist_range));
-            problem.SetParameterLowerBound(&distance_mm, 0, distance_mm * (1.0 + dist_range));
+            problem.SetParameterUpperBound(&distance_mm, 0, distance_mm * (1.0 + dist_range));
         }
 
         if (!data.refine_beam_center) {