From cd561dd3d989946b84a4a089de00957d4bb9eda9 Mon Sep 17 00:00:00 2001
From: Filip Leonarski <filip.leonarski@psi.ch>
Date: Wed, 4 Mar 2026 09:00:43 +0100
Subject: [PATCH] ScaleAndMerge: rotation wedge can be refined

---
 image_analysis/scale_merge/ScaleAndMerge.cpp | 19 +++++++++++++------
 image_analysis/scale_merge/ScaleAndMerge.h   |  2 ++
 2 files changed, 15 insertions(+), 6 deletions(-)

diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp
index a45264bc..2d813178 100644
--- a/image_analysis/scale_merge/ScaleAndMerge.cpp
+++ b/image_analysis/scale_merge/ScaleAndMerge.cpp
@@ -121,7 +121,6 @@ namespace {
               weight_(SafeInv(sigma_obs, 1.0)),
               delta_phi_(r.delta_phi_deg),
               lp_(SafeInv(r.rlp, 1.0)),
-              half_wedge_(wedge_deg / 2.0),
               c1_(r.zeta / std::sqrt(2.0)),
               partiality_(r.partiality) {
         }
@@ -130,11 +129,13 @@ namespace {
         bool operator()(const T *const G,
                         const T *const mosaicity,
                         const T *const Itrue,
+                        const T *const wedge,
                         T *residual) const {
             T partiality;
             if (mosaicity[0] >= 0.0) {
-                const T arg_plus = T(delta_phi_ + half_wedge_) * T(c1_) / mosaicity[0];
-                const T arg_minus = T(delta_phi_ - half_wedge_) * T(c1_) / mosaicity[0];
+                const T half_wedge = wedge[0] / T(2.0);
+                const T arg_plus = T(delta_phi_ + half_wedge) * T(c1_) / mosaicity[0];
+                const T arg_minus = T(delta_phi_ - half_wedge) * T(c1_) / mosaicity[0];
                 partiality = (ceres::erf(arg_plus) - ceres::erf(arg_minus)) / T(2.0);
             } else
                 partiality = T(1.0);
@@ -148,7 +149,6 @@ namespace {
         double weight_;
         double delta_phi_;
         double lp_;
-        double half_wedge_;
         double c1_;
         double partiality_;
     };
@@ -275,17 +275,20 @@ namespace {
             }
         }
 
+        double wedge = opt.wedge_deg.value_or(0.0);
+
         std::vector<bool> is_valid_hkl_slot(nhkl, false);
         for (const auto &o: obs) {
             switch (opt.partiality_model) {
                 case ScaleMergeOptions::PartialityModel::Rotation: {
-                    auto *cost = new ceres::AutoDiffCostFunction<IntensityRotResidual, 1, 1, 1, 1>(
+                    auto *cost = new ceres::AutoDiffCostFunction<IntensityRotResidual, 1, 1, 1, 1, 1>(
                         new IntensityRotResidual(*o.r, o.sigma, opt.wedge_deg.value_or(0.0)));
                     problem.AddResidualBlock(cost,
                                              nullptr,
                                              &g[o.img_id],
                                              &mosaicity[o.img_id],
-                                             &Itrue[o.hkl_slot]);
+                                             &Itrue[o.hkl_slot],
+                                             &wedge);
                 }
                 break;
                 case ScaleMergeOptions::PartialityModel::Still: {
@@ -374,6 +377,10 @@ namespace {
                     problem.SetParameterUpperBound(&mosaicity[i], 0, opt.mosaicity_max_deg);
                 }
             }
+            if (!opt.refine_wedge)
+                problem.SetParameterBlockConstant(&wedge);
+            else
+                problem.SetParameterLowerBound(&wedge, 0, 0.0);
         }
 
         // use all available threads
diff --git a/image_analysis/scale_merge/ScaleAndMerge.h b/image_analysis/scale_merge/ScaleAndMerge.h
index 9710ce05..d154e941 100644
--- a/image_analysis/scale_merge/ScaleAndMerge.h
+++ b/image_analysis/scale_merge/ScaleAndMerge.h
@@ -48,6 +48,8 @@ struct ScaleMergeOptions {
 
     int64_t image_cluster = 1;
 
+    bool refine_wedge = false;
+
     enum class PartialityModel {Fixed, Rotation, Unity, Still} partiality_model = PartialityModel::Fixed;
 };