From 068ea2eea768c9fbf230c45576bf0b8b3285b786 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 8 Feb 2026 14:36:06 +0100 Subject: [PATCH] ScaleAndMerge: Further improvement --- image_analysis/scale_merge/ScaleAndMerge.cpp | 90 ++++++++++++++++---- image_analysis/scale_merge/ScaleAndMerge.h | 22 +++-- 2 files changed, 90 insertions(+), 22 deletions(-) diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index 412d7b8b..bc85d816 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -126,7 +126,8 @@ struct IntensityResidual { half_wedge_rad_(wedge_deg * M_PI / 180.0 / 2.0), c1_(std::sqrt(2.0) * SafeInv(static_cast(r.zeta), 1.0)), s2_(s2), - refine_b_(refine_b) {} + refine_b_(refine_b), + partiality_model_(partiality_model) {} template bool operator()(const T* const log_k, @@ -138,7 +139,7 @@ struct IntensityResidual { const T Itrue = ceres::exp(log_Itrue[0]); T partiality = T(1.0); - if (partiality_model) { + if (partiality_model_) { // partiality = 0.5 * [ erf( (Δφ+Δω/2) * (sqrt(2)*η/ζ) ) - erf( (Δφ-Δω/2) * (sqrt(2)*η/ζ) ) ] const T arg_plus = T(delta_phi_rad_ + half_wedge_rad_) * (T(c1_) * mosaicity_rad[0]); const T arg_minus = T(delta_phi_rad_ - half_wedge_rad_) * (T(c1_) * mosaicity_rad[0]); @@ -162,7 +163,21 @@ struct IntensityResidual { double c1_; double s2_; bool refine_b_; - bool partiality_model; + bool partiality_model_; +}; + +struct ScaleRegularizationResidual { + explicit ScaleRegularizationResidual(double sigma_k) + : inv_sigma_(SafeInv(sigma_k, 1.0)) {} + + template + bool operator()(const T* const log_k, T* residual) const { + const T k = ceres::exp(log_k[0]); + residual[0] = (k - T(1.0)) * T(inv_sigma_); + return true; + } + + double inv_sigma_; }; } // namespace @@ -257,8 +272,17 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector& ob std::vector b(nimg, 0.0); std::vector log_Itrue(nhkl, 0.0); - // Global mosaicity (radians) - double mosaicity_rad = opt.mosaicity_init_rad; + auto deg2rad = [](double deg) { return deg * M_PI / 180.0; }; + + // Mosaicity: either global or per-image (always stored in radians internally) + std::vector mosaicity_rad; + double mosaicity_rad_global = deg2rad(opt.mosaicity_init_deg); + + if (opt.mosaicity_per_image) { + mosaicity_rad.assign(nimg, deg2rad(opt.mosaicity_init_deg)); + } else { + mosaicity_rad_global = deg2rad(opt.mosaicity_init_deg); + } // Initialize Itrue from per-HKL median of observed intensities (rough; model contains partiality/LP) { @@ -282,26 +306,36 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector& ob ceres::Problem problem; - std::unique_ptr loss; + std::unique_ptr loss; if (opt.use_huber_loss) { loss = std::make_unique(opt.huber_delta); } - bool partiality_model = opt.wedge_deg > 0.0; + const bool partiality_model = opt.wedge_deg > 0.0; for (const auto& o : obs) { + double* mosaicity_block = opt.mosaicity_per_image ? &mosaicity_rad[o.img_slot] : &mosaicity_rad_global; + auto* cost = new ceres::AutoDiffCostFunction( - new IntensityResidual(*o.r, o.sigma, opt.wedge_deg, o.s2, opt.refine_b_factor, - partiality_model)); + new IntensityResidual(*o.r, o.sigma, opt.wedge_deg, o.s2, opt.refine_b_factor, partiality_model)); problem.AddResidualBlock(cost, loss.get(), &log_k[o.img_slot], &b[o.img_slot], - &mosaicity_rad, + mosaicity_block, &log_Itrue[o.hkl_slot]); } + // Optional Kabsch-like regularization for k: (k - 1)/sigma + if (opt.regularize_scale_to_one) { + for (int i = 0; i < nimg; ++i) { + auto* rcost = new ceres::AutoDiffCostFunction( + new ScaleRegularizationResidual(opt.scale_regularization_sigma)); + problem.AddResidualBlock(rcost, nullptr, &log_k[i]); + } + } + // Fix gauge freedom: anchor first image scale to 1.0 if (opt.fix_first_image_scale && nimg > 0) { log_k[0] = 0.0; @@ -322,13 +356,27 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector& ob } } + // Mosaicity refinement + bounds (degrees -> radians) if (!opt.refine_mosaicity) { - problem.SetParameterBlockConstant(&mosaicity_rad); + if (opt.mosaicity_per_image) { + for (int i = 0; i < nimg; ++i) + problem.SetParameterBlockConstant(&mosaicity_rad[i]); + } else { + problem.SetParameterBlockConstant(&mosaicity_rad_global); + } } else { - if (opt.mosaicity_min_rad) - problem.SetParameterLowerBound(&mosaicity_rad, 0, *opt.mosaicity_min_rad); - if (opt.mosaicity_max_rad) - problem.SetParameterUpperBound(&mosaicity_rad, 0, *opt.mosaicity_max_rad); + const std::optional min_rad = opt.mosaicity_min_deg ? std::optional(deg2rad(*opt.mosaicity_min_deg)) : std::nullopt; + const std::optional max_rad = opt.mosaicity_max_deg ? std::optional(deg2rad(*opt.mosaicity_max_deg)) : std::nullopt; + + if (opt.mosaicity_per_image) { + for (int i = 0; i < nimg; ++i) { + if (min_rad) problem.SetParameterLowerBound(&mosaicity_rad[i], 0, *min_rad); + if (max_rad) problem.SetParameterUpperBound(&mosaicity_rad[i], 0, *max_rad); + } + } else { + if (min_rad) problem.SetParameterLowerBound(&mosaicity_rad_global, 0, *min_rad); + if (max_rad) problem.SetParameterUpperBound(&mosaicity_rad_global, 0, *max_rad); + } } ceres::Solver::Options options; @@ -346,6 +394,15 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector& ob out.image_b_factor[i] = opt.refine_b_factor ? b[i] : 0.0; } + // Export mosaicity (degrees) to result + if (opt.mosaicity_per_image) { + out.mosaicity_deg.resize(nimg); + for (int i = 0; i < nimg; ++i) + out.mosaicity_deg[i] = mosaicity_rad[i] * 180.0 / M_PI; + } else { + out.mosaicity_deg = { mosaicity_rad_global * 180.0 / M_PI }; + } + // Reverse maps for merged output std::vector slotToHKL(nhkl); for (const auto& kv : hklToSlot) { @@ -367,7 +424,8 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector& ob const double Itrue = std::exp(log_Itrue[h]); const double zeta = static_cast(o.r->zeta); - const double c1 = std::sqrt(2.0) * (mosaicity_rad / zeta); + const double mosa_i = opt.mosaicity_per_image ? mosaicity_rad[i] : mosaicity_rad_global; + const double c1 = std::sqrt(2.0) * (mosa_i / zeta); const double delta_phi_rad = static_cast(o.r->delta_phi) * M_PI / 180.0; const double partiality = partiality_model ? diff --git a/image_analysis/scale_merge/ScaleAndMerge.h b/image_analysis/scale_merge/ScaleAndMerge.h index 93806712..848c7509 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.h +++ b/image_analysis/scale_merge/ScaleAndMerge.h @@ -35,16 +35,23 @@ struct ScaleMergeOptions { // If false, keep them separate by including a sign flag in the HKL key. bool merge_friedel = true; - // --- New: parameters for Kabsch(XDS)-style partiality model --- + // --- Kabsch(XDS)-style partiality model --- // Rotation range (wedge) used in partiality calculation. + // Set to 0 to disable partiality correction. double wedge_deg = 1.0; - // Refine global mosaicity (sigma, in radians). If false, mosaicity is fixed. + // --- Mosaicity (user input in degrees; internally converted to radians) --- + // If true, refine mosaicity. If mosaicity_per_image==true, refine one value per image. bool refine_mosaicity = true; - double mosaicity_init_rad = 0.003; // ~0.17 deg - std::optional mosaicity_min_rad = 1e-5; - std::optional mosaicity_max_rad = 0.2; // generous upper bound + bool mosaicity_per_image = true; + double mosaicity_init_deg = 0.17; // ~0.003 rad + std::optional mosaicity_min_deg = 1e-3; + std::optional mosaicity_max_deg = 20.0; + + // --- Optional: regularize per-image scale k towards 1 (Kabsch-like) --- + bool regularize_scale_to_one = false; + double scale_regularization_sigma = 0.05; }; struct MergedReflection { @@ -60,7 +67,10 @@ struct ScaleMergeResult { std::vector image_scale_k; std::vector image_b_factor; std::vector image_ids; + + // If mosaicity_per_image==true, one value per image (degrees). Otherwise size is 1 (global). + std::vector mosaicity_deg; }; ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector& observations, - const ScaleMergeOptions& opt = {}); + const ScaleMergeOptions& opt = {}); \ No newline at end of file