XtalOptimizer: Store 1/lambda (not lambda) + add check for very small lambda + set beam center as one parameter block + fix poni rotation convention (rot2 was likely wrong)

image_analysis/geom_refinement/XtalOptimizer.cpp

@@ -16,18 +16,20 @@ struct XtalResidual {
                  double exp_l,
                  gemmi::CrystalSystem symmetry)
         : obs_x(x), obs_y(y),
-          lambda(lambda),
+          inv_lambda(1.0/lambda),
           pixel_size(pixel_size),
           exp_h(exp_h),
           exp_k(exp_k),
           exp_l(exp_l),
           angle_rad(angle_rad),
           symmetry(symmetry) {
+        if (std::fabs(lambda) < 1e-6)
+            throw JFJochException(JFJochExceptionCategory::InputParameterInvalid,
+                "Lambda cannot be close to zero");
     }
 
     template<typename T>
-    bool operator()(const T *const beam_x,
-                    const T *const beam_y,
+    bool operator()(const T *const beam,
                     const T *const distance_mm,
                     const T *const detector_rot,
                     const T *const rotation_axis,
@@ -47,12 +49,12 @@ struct XtalResidual {
         const T s1 = ceres::sin(rot1);
 
         // Rx(-rot2): rotation around X-axis with inverted sign (PyFAI left-handed)
-        const T c2 = ceres::cos(-rot2);
-        const T s2 = ceres::sin(-rot2);
+        const T c2 = ceres::cos(rot2);
+        const T s2 = ceres::sin(rot2);
 
         // Detector coordinates in mm
-        const T det_x = (T(obs_x) - beam_x[0]) * T(pixel_size);
-        const T det_y = (T(obs_y) - beam_y[0]) * T(pixel_size);
+        const T det_x = (T(obs_x) - beam[0]) * T(pixel_size);
+        const T det_y = (T(obs_y) - beam[1]) * T(pixel_size);
         const T det_z = T(distance_mm[0]);
 
         // Apply Ry(rot1) first: rotate around Y
@@ -62,18 +64,17 @@ struct XtalResidual {
 
         // Then apply Rx(-rot2): rotate around X
         const T x = t1_x;
-        const T y = c2 * t1_y - s2 * t1_z;
-        const T z = s2 * t1_y + c2 * t1_z;
+        const T y = c2 * t1_y + s2 * t1_z;
+        const T z = -s2 * t1_y + c2 * t1_z;
 
         // convert to recip space
         const T lab_norm = ceres::sqrt(x * x + y * y + z * z);
         const T inv_norm = T(1) / lab_norm;
-        const T inv_lambda = T(1) / T(lambda);
 
         T recip_raw[3];
-        recip_raw[0] = x * inv_norm * inv_lambda;
-        recip_raw[1] = y * inv_norm * inv_lambda;
-        recip_raw[2] = (z * inv_norm - T(1.0)) * inv_lambda;
+        recip_raw[0] = x * inv_norm * T(inv_lambda);
+        recip_raw[1] = y * inv_norm * T(inv_lambda);
+        recip_raw[2] = (z * inv_norm - T(1.0)) * T(inv_lambda);
 
         // Apply goniometer "back-to-start" rotation:
         // brings observed reciprocal from image orientation into reference crystal frame
@@ -178,7 +179,7 @@ struct XtalResidual {
     }
 
     const double obs_x, obs_y;
-    const double lambda;
+    const double inv_lambda;
     const double pixel_size;
     const double exp_h;
     const double exp_k;
@@ -510,17 +511,16 @@ bool XtalOptimizerInternal(XtalOptimizerData &data,
         double beta = data.latt.GetUnitCell().beta;
 
         // Initial guess for the parameters
-        double beam_x = data.geom.GetBeamX_pxl();
-        double beam_y = data.geom.GetBeamY_pxl();
+        double beam[2] = {data.geom.GetBeamX_pxl(), data.geom.GetBeamY_pxl()};
         double distance_mm = data.geom.GetDetectorDistance_mm();
 
         double detector_rot[2] = {data.geom.GetPoniRot1_rad(), data.geom.GetPoniRot2_rad()};
 
         ceres::Problem problem;
 
-double latt_vec0[3] = {0.0, 0.0, 0.0};
-double latt_vec1[3] = {0.0, 0.0, 0.0};
-double latt_vec2[3] = {0.0, 0.0, 0.0};
+        double latt_vec0[3] = {0.0, 0.0, 0.0};
+        double latt_vec1[3] = {0.0, 0.0, 0.0};
+        double latt_vec2[3] = {0.0, 0.0, 0.0};
 
         double rot_vec[3] = {1, 0, 0};
 
@@ -591,7 +591,7 @@ double latt_vec2[3] = {0.0, 0.0, 0.0};
                 continue;
 
             problem.AddResidualBlock(
-                new ceres::AutoDiffCostFunction<XtalResidual, 3, 1, 1, 1, 2, 3, 3, 3, 3>(
+                new ceres::AutoDiffCostFunction<XtalResidual, 3, 2, 1, 2, 3, 3, 3, 3>(
                     new XtalResidual(pt.x, pt.y,
                                      data.geom.GetWavelength_A(),
                                      data.geom.GetPixelSize_mm(),
@@ -599,8 +599,7 @@ double latt_vec2[3] = {0.0, 0.0, 0.0};
                                      h, k, l,
                                      data.crystal_system)),
                 nullptr,
-                &beam_x,
-                &beam_y,
+                beam,
                 &distance_mm,
                 detector_rot,
                 rot_vec,
@@ -621,10 +620,8 @@ double latt_vec2[3] = {0.0, 0.0, 0.0};
             problem.SetParameterUpperBound(&distance_mm, 0, distance_mm * (1.0 + dist_range));
         }
 
-        if (!data.refine_beam_center) {
-            problem.SetParameterBlockConstant(&beam_x);
-            problem.SetParameterBlockConstant(&beam_y);
-        }
+        if (!data.refine_beam_center)
+            problem.SetParameterBlockConstant(beam);
 
         if (!data.refine_detector_angles) {
             problem.SetParameterBlockConstant(detector_rot);
@@ -684,9 +681,9 @@ double latt_vec2[3] = {0.0, 0.0, 0.0};
         ceres::Solve(options, &problem, &summary);
 
         if (data.refine_beam_center) {
-            data.beam_corr_x = data.geom.GetBeamX_pxl() - beam_x;
-            data.beam_corr_y = data.geom.GetBeamY_pxl() - beam_y;
-            data.geom.BeamX_pxl(beam_x).BeamY_pxl(beam_y);
+            data.beam_corr_x = data.geom.GetBeamX_pxl() - beam[0];
+            data.beam_corr_y = data.geom.GetBeamY_pxl() - beam[1];
+            data.geom.BeamX_pxl(beam[0]).BeamY_pxl(beam[1]);
         }
 
         if (data.refine_distance_mm)