From 89d827a24f70f9cf077ebf77d824bbc71036f305 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 19 Feb 2026 10:29:38 +0100 Subject: [PATCH] jfjoch_process: Add resolution limit --- image_analysis/scale_merge/ScaleAndMerge.cpp | 7 +++++++ image_analysis/scale_merge/ScaleAndMerge.h | 1 + tools/jfjoch_process.cpp | 14 ++++++++++---- 3 files changed, 18 insertions(+), 4 deletions(-) diff --git a/image_analysis/scale_merge/ScaleAndMerge.cpp b/image_analysis/scale_merge/ScaleAndMerge.cpp index 8c6c63ca..30733018 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.cpp +++ b/image_analysis/scale_merge/ScaleAndMerge.cpp @@ -268,6 +268,9 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector &ob if (!std::isfinite(r.I)) continue; + if (opt.d_min_limit_A > 0.0 && d < opt.d_min_limit_A) + continue; + if (!std::isfinite(r.zeta) || r.zeta <= 0.0f) continue; if (!std::isfinite(r.rlp) || r.rlp == 0.0f) @@ -582,6 +585,8 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector &ob for (int h = 0; h < nhkl; ++h) { const auto d = static_cast(out.merged[h].d); if (std::isfinite(d) && d > 0.0f) { + if (opt.d_min_limit_A > 0.0 && d < static_cast(opt.d_min_limit_A)) + continue; stat_d_min = std::min(stat_d_min, d); stat_d_max = std::max(stat_d_max, d); } @@ -600,6 +605,8 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector &ob for (int h = 0; h < nhkl; ++h) { const auto d = static_cast(out.merged[h].d); if (std::isfinite(d) && d > 0.0f) { + if (opt.d_min_limit_A > 0.0 && d < static_cast(opt.d_min_limit_A)) + continue; auto s = stat_shells.GetShell(d); if (s.has_value()) hkl_shell[h] = s.value(); diff --git a/image_analysis/scale_merge/ScaleAndMerge.h b/image_analysis/scale_merge/ScaleAndMerge.h index 7cf9a3c4..439ec25d 100644 --- a/image_analysis/scale_merge/ScaleAndMerge.h +++ b/image_analysis/scale_merge/ScaleAndMerge.h @@ -48,6 +48,7 @@ struct ScaleMergeOptions { bool smoothen_g = true; bool smoothen_mos = true; + double d_min_limit_A = 0.0; }; struct MergedReflection { diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 646392f0..208a5723 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -39,6 +39,7 @@ void print_usage(Logger &logger) { logger.Info(" -F Use FFT indexing algorithm (default: Auto)"); logger.Info(" -x No least-square beam center refinement"); logger.Info(" -d High resolution limit for spot finding (default: 1.5)"); + 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"); @@ -67,6 +68,8 @@ 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; @@ -79,7 +82,7 @@ int main(int argc, char **argv) { } int opt; - while ((opt = getopt(argc, argv, "o:N:s:e:vR::Fxd:S:MLm:A")) != -1) { + while ((opt = getopt(argc, argv, "o:N:s:e:vR::Fxd:S:MLm:AD:")) != -1) { switch (opt) { case 'o': output_prefix = optarg; @@ -96,6 +99,9 @@ int main(int argc, char **argv) { case 'v': verbose = true; break; + case 'd': + resolution_limit = atof(optarg); + break; case 'R': rotation_indexing = true; if (optarg) rotation_indexing_range = atof(optarg); @@ -106,8 +112,9 @@ int main(int argc, char **argv) { case 'x': refine_beam_center = false; break; - case 'd': + case 'D': d_high = atof(optarg); + logger.Info("High resolution limit for scaling/merging set to {:.2f} A", d_high); break; case 'S': space_group_number = atoi(optarg); @@ -425,6 +432,7 @@ int main(int argc, char **argv) { 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; switch (mosaicity_refinement_mode) { case MosaicityRefinementMode::None: @@ -450,8 +458,6 @@ int main(int argc, char **argv) { double scale_time = std::chrono::duration(scale_end - scale_start).count(); if (scale_result) { - -// ... existing code ... logger.Info("Scaling completed in {:.2f} s ({} unique reflections, {} images)", scale_time, scale_result->merged.size(),