From 43c340219892901ec33c4bb00d196d30605050b8 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Fri, 20 Feb 2026 09:28:50 +0100 Subject: [PATCH] ScaleMerge: Simplify by removing log residual --- image_analysis/scale_merge/ScaleAndMerge.cpp | 95 ++------------------ image_analysis/scale_merge/ScaleAndMerge.h | 1 - tools/jfjoch_process.cpp | 25 +++--- 3 files changed, 18 insertions(+), 103 deletions(-) diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index 30733018..9aeab336 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -115,59 +115,6 @@ namespace { return key; } - /// CrystFEL-like log-scaling residual - /// - /// residual = w * [ ln(I_obs) - ln(G) - ln(partiality) - ln(lp) - ln(I_true) ] - /// - /// Only observations with I_obs > 0 should be fed in (the caller skips the rest). - /// G and I_true are constrained to be positive via Ceres lower bounds. - struct LogIntensityResidual { - LogIntensityResidual(const Reflection &r, double sigma_obs, double wedge_deg, bool refine_partiality) - : log_Iobs_(std::log(std::max(static_cast(r.I), 1e-30))), - weight_(SafeInv(sigma_obs / std::max(static_cast(r.I), 1e-30), 1.0)), - delta_phi_(r.delta_phi_deg), - log_lp_(std::log(std::max(SafeInv(static_cast(r.rlp), 1.0), 1e-30))), - half_wedge_(wedge_deg / 2.0), - c1_(r.zeta / std::sqrt(2.0)), - partiality_(r.partiality), - refine_partiality_(refine_partiality) { - } - - template - bool operator()(const T *const G, - const T *const mosaicity, - const T *const Itrue, - T *residual) const { - T partiality; - if (refine_partiality_ && 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]; - partiality = (ceres::erf(arg_plus) - ceres::erf(arg_minus)) / T(2.0); - } else { - partiality = T(partiality_); - } - - // Clamp partiality away from zero so log is safe - const T min_p = T(1e-30); - if (partiality < min_p) - partiality = min_p; - - // ln(I_pred) = ln(G) + ln(partiality) + ln(lp) + ln(Itrue) - const T log_Ipred = ceres::log(G[0]) + ceres::log(partiality) + T(log_lp_) + ceres::log(Itrue[0]); - residual[0] = (log_Ipred - T(log_Iobs_)) * T(weight_); - return true; - } - - double log_Iobs_; - double weight_; // w_i ≈ I_obs / sigma_obs (relative weight in log-space) - double delta_phi_; - double log_lp_; - double half_wedge_; - double c1_; - double partiality_; - bool refine_partiality_; - }; - struct IntensityResidual { IntensityResidual(const Reflection &r, double sigma_obs, double wedge_deg, bool refine_partiality) : Iobs_(static_cast(r.I)), @@ -362,32 +309,15 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector &ob for (const auto &o: obs) { size_t mos_slot = opt.per_image_mosaicity ? o.img_slot : 0; + auto *cost = new ceres::AutoDiffCostFunction( + new IntensityResidual(*o.r, o.sigma, opt.wedge_deg, refine_partiality)); - if (opt.log_scaling_residual) { - // Log residual requires positive I_obs - if (o.r->I <= 0.0f) - continue; - - auto *cost = new ceres::AutoDiffCostFunction( - new LogIntensityResidual(*o.r, o.sigma, opt.wedge_deg, refine_partiality)); - - problem.AddResidualBlock(cost, - nullptr, - &g[o.img_slot], - &mosaicity[mos_slot], - &Itrue[o.hkl_slot]); - is_valid_hkl_slot[o.hkl_slot] = true; - } else { - auto *cost = new ceres::AutoDiffCostFunction( - new IntensityResidual(*o.r, o.sigma, opt.wedge_deg, refine_partiality)); - - problem.AddResidualBlock(cost, - nullptr, - &g[o.img_slot], - &mosaicity[mos_slot], - &Itrue[o.hkl_slot]); - is_valid_hkl_slot[o.hkl_slot] = true; - } + problem.AddResidualBlock(cost, + nullptr, + &g[o.img_slot], + &mosaicity[mos_slot], + &Itrue[o.hkl_slot]); + is_valid_hkl_slot[o.hkl_slot] = true; } if (opt.smoothen_g) { @@ -418,15 +348,6 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector &ob for (int i = 0; i < nimg; ++i) problem.SetParameterLowerBound(&g[i], 0, 1e-12); - - // For log residual, Itrue must stay positive - if (opt.log_scaling_residual) { - for (int h = 0; h < nhkl; ++h) { - if (is_valid_hkl_slot[h]) - problem.SetParameterLowerBound(&Itrue[h], 0, 1e-12); - } - } - // Mosaicity refinement + bounds if (!opt.refine_mosaicity) { for (int i = 0; i < (opt.per_image_mosaicity ? nimg : 1); ++i) diff --git a/image_analysis/scale_merge/ScaleAndMerge.h b/image_analysis/scale_merge/ScaleAndMerge.h index 439ec25d..8d766785 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.h +++ b/image_analysis/scale_merge/ScaleAndMerge.h @@ -42,7 +42,6 @@ struct ScaleMergeOptions { // --- Optional: regularize per-image scale k towards 1 (Kabsch-like) --- bool regularize_scale_to_one = true; double scale_regularization_sigma = 0.05; - bool log_scaling_residual = false; double min_partiality_for_merge = 0.2; bool smoothen_g = true; diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 208a5723..b24b6af2 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -42,7 +42,6 @@ void print_usage(Logger &logger) { logger.Info(" -D High resolution limit for scaling/merging (default: 0.0; no limit)"); logger.Info(" -S Space group number"); logger.Info(" -M Scale and merge (refine mosaicity) and write scaled.hkl + image.dat"); - logger.Info(" -L Use log-scaling residual"); logger.Info(" -m Mosaicity refinement none|fixed|image (default: image)"); logger.Info(" -A Anomalous mode (don't merge Friedel pairs)"); } @@ -68,13 +67,11 @@ int main(int argc, char **argv) { bool anomalous_mode = false; std::optional space_group_number; - double resolution_limit = 0.0; - - bool log_residual = false; enum class MosaicityRefinementMode { None, Fixed, Image }; MosaicityRefinementMode mosaicity_refinement_mode = MosaicityRefinementMode::Image; - float d_high = 1.5; + float d_min_spot_finding = 1.5; + std::optional d_min_scale_merge; if (argc == 1) { print_usage(logger); @@ -82,7 +79,7 @@ int main(int argc, char **argv) { } int opt; - while ((opt = getopt(argc, argv, "o:N:s:e:vR::Fxd:S:MLm:AD:")) != -1) { + while ((opt = getopt(argc, argv, "o:N:s:e:vR::Fxd:S:Mm:AD:")) != -1) { switch (opt) { case 'o': output_prefix = optarg; @@ -100,7 +97,7 @@ int main(int argc, char **argv) { verbose = true; break; case 'd': - resolution_limit = atof(optarg); + d_min_spot_finding = atof(optarg); break; case 'R': rotation_indexing = true; @@ -113,8 +110,8 @@ int main(int argc, char **argv) { refine_beam_center = false; break; case 'D': - d_high = atof(optarg); - logger.Info("High resolution limit for scaling/merging set to {:.2f} A", d_high); + d_min_scale_merge = atof(optarg); + logger.Info("High resolution limit for scaling/merging set to {:.2f} A", d_min_spot_finding); break; case 'S': space_group_number = atoi(optarg); @@ -122,9 +119,6 @@ int main(int argc, char **argv) { case 'M': run_scaling = true; break; - case 'L': - log_residual = true; - break; case 'A': anomalous_mode = true; break; @@ -213,7 +207,9 @@ int main(int argc, char **argv) { SpotFindingSettings spot_settings; spot_settings.enable = true; spot_settings.indexing = true; - spot_settings.high_resolution_limit = d_high; + spot_settings.high_resolution_limit = d_min_spot_finding; + if (d_min_scale_merge > 0) + spot_settings.high_resolution_limit = d_min_spot_finding; // Initialize Analysis Components PixelMask pixel_mask = dataset->pixel_mask; @@ -430,9 +426,8 @@ int main(int argc, char **argv) { scale_opts.refine_mosaicity = true; scale_opts.max_num_iterations = 500; scale_opts.max_solver_time_s = 240.0; // generous cutoff for now - scale_opts.log_scaling_residual = log_residual; scale_opts.merge_friedel = !anomalous_mode; - scale_opts.d_min_limit_A = d_high; + scale_opts.d_min_limit_A = d_min_scale_merge.value_or(0.0); switch (mosaicity_refinement_mode) { case MosaicityRefinementMode::None: