jfjoch_process: Add resolution limit

image_analysis/scale_merge/ScaleAndMerge.cpp

@@ -268,6 +268,9 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<Reflection> &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<Reflection> &ob
         for (int h = 0; h < nhkl; ++h) {
             const auto d = static_cast<float>(out.merged[h].d);
             if (std::isfinite(d) && d > 0.0f) {
+                if (opt.d_min_limit_A > 0.0 && d < static_cast<float>(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<Reflection> &ob
             for (int h = 0; h < nhkl; ++h) {
                 const auto d = static_cast<float>(out.merged[h].d);
                 if (std::isfinite(d) && d > 0.0f) {
+                    if (opt.d_min_limit_A > 0.0 && d < static_cast<float>(opt.d_min_limit_A))
+                        continue;
                     auto s = stat_shells.GetShell(d);
                     if (s.has_value())
                         hkl_shell[h] = s.value();

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 {

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<num>    High resolution limit for spot finding (default: 1.5)");
+    logger.Info("   -D<num>    High resolution limit for scaling/merging (default: 0.0; no limit)");
     logger.Info("   -S<num>    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<int> 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<double>(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(),