jfjoch_process: reflection selection for scaling

image_analysis/scale_merge/ScaleAndMerge.cpp

@@ -235,6 +235,7 @@ namespace {
         int img_id = 0;
         int hkl_slot = -1;
         double sigma = 0.0;
+        mutable bool selected = true;
     };
 
     struct CorrectedObs {
@@ -243,6 +244,104 @@ namespace {
         double sigma_corr;
     };
 
+    void select_reflections_by_quasi_random(const std::vector<ObsRef> &obs,
+        const ScaleMergeOptions &opt,
+        std::vector<bool> &hkl_selected,
+        double Isigma_cutoff = 1.0,
+        int min_per_bin = 10000,
+        int max_per_bin = 80000,
+        int n_resolution_bins = 20) {
+        float stat_d_min = std::numeric_limits<float>::max();
+        float stat_d_max = 0.0f;
+
+        struct HKLStats {
+            int n = 0;
+            float d = std::numeric_limits<float>::max();
+            int shell_id = 0;
+        };
+        const int nhkl = static_cast<int>(hkl_selected.size());
+        std::vector<HKLStats> per_hkl(nhkl);
+        int reflection_above_cutoff = 0;
+
+        for (const auto &o: obs) {
+            if (o.r->I / o.r->sigma < Isigma_cutoff) {
+                o.selected = false;
+                continue;
+            }
+
+            const auto d = o.r->d;
+            reflection_above_cutoff += 1;
+            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);
+                auto &hs = per_hkl[o.hkl_slot];
+                hs.n += 1;
+                hs.d = d;
+            }
+        }
+
+        std::cout << "Reflections of I/sigma > " << Isigma_cutoff << " : " << reflection_above_cutoff << std::endl;
+        if (reflection_above_cutoff < min_per_bin * n_resolution_bins) {
+            std::cout << "No additional selection applied before scaling" << std::endl;
+            return;
+        }
+
+        if (stat_d_min < stat_d_max && stat_d_min > 0.0f) {
+            const float d_min_pad = stat_d_min * 0.999f;
+            const float d_max_pad = stat_d_max * 1.001f;
+            ResolutionShells scaling_shells(d_min_pad, d_max_pad, n_resolution_bins);
+
+            for (int h = 0; h < nhkl; ++h) {
+                const auto d = per_hkl[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 = scaling_shells.GetShell(d);
+                    if (s.has_value())
+                        per_hkl[h].shell_id = s.value();
+                }
+            }
+
+            // Accumulators per shell
+            struct ShellAccum {
+                int obs_unique = 0;
+                int obs_total = 0;
+                bool selected = true;
+            };
+            std::vector<ShellAccum> shell_acc(n_resolution_bins);
+
+            for (int h = 0; h < nhkl; ++h) {
+                auto &sa = shell_acc[per_hkl[h].shell_id];
+                if (sa.obs_unique > min_per_bin * 1.2 || sa.obs_total > max_per_bin)
+                    hkl_selected[h] = false;
+                else
+                    sa.obs_unique += 1;
+                    sa.obs_total += per_hkl[h].n;
+            }
+
+            const auto shell_min_res = scaling_shells.GetShellMinRes();
+
+            std::cout << "| d-mean | n_refl_tot | n_refl_uni | selection |" << std::endl;
+            for (int n=0; n < n_resolution_bins; ++n) {
+                if (shell_acc[n].obs_unique < min_per_bin)
+                    shell_acc[n].selected = false;
+                if (shell_min_res[n] <= 0.0f) continue;
+                std::cout << std::setw(8) << std::fixed << std::setprecision(3) << shell_min_res[n]
+                          << std::setw(12) << std::fixed << std::setprecision(0) << shell_acc[n].obs_unique
+                          << std::setw(12) << shell_acc[n].obs_total
+                          << "      " << shell_acc[n].selected << std::endl;
+            }
+
+            for (int h = 0; h < nhkl; ++h) {
+                auto &sa = shell_acc[per_hkl[h].shell_id];
+                if (!sa.selected)
+                    hkl_selected[h] = false;
+            }
+       }
+    }
+
     void scale(const ScaleMergeOptions &opt,
                std::vector<double> &g,
                std::vector<double> &mosaicity,
@@ -250,10 +349,12 @@ namespace {
                const std::vector<uint8_t> &image_slot_used,
                bool rotation_crystallography,
                size_t nhkl,
-               const std::vector<ObsRef> &obs) {
+               const std::vector<ObsRef> &obs,
+               bool selection) {
         ceres::Problem problem;
 
         std::vector<double> Itrue(nhkl, 0.0);
+        std::vector<bool> hkl_selected(nhkl, true);
 
         // Initialize Itrue from per-HKL median of observed intensities
         {
@@ -277,8 +378,12 @@ namespace {
 
         double wedge = opt.wedge_deg.value_or(0.0);
 
+        select_reflections_by_quasi_random(obs, opt, hkl_selected);
+
         std::vector<bool> is_valid_hkl_slot(nhkl, false);
         for (const auto &o: obs) {
+            if (!o.selected) continue;
+            if (!hkl_selected[o.hkl_slot]) continue;
             switch (opt.partiality_model) {
                 case ScaleMergeOptions::PartialityModel::Rotation: {
                     auto *cost = new ceres::AutoDiffCostFunction<IntensityRotResidual, 1, 1, 1, 1, 1>(
@@ -398,6 +503,7 @@ namespace {
         options.function_tolerance = 1e-4;
 
         ceres::Solver::Summary summary;
+        std::cout << "Now start the ceres-solver with residual blocks: " << problem.NumResidualBlocks() << std::endl;
         ceres::Solve(options, &problem, &summary);
         std::cout << summary.FullReport() << std::endl;
     }
@@ -694,6 +800,7 @@ namespace {
                 o.img_id = img_id;
                 o.hkl_slot = hkl_slot;
                 o.sigma = sigma;
+                o.selected = true;
                 obs.push_back(o);
             }
         }
@@ -707,6 +814,7 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<std::vector<Ref
         throw std::invalid_argument("image_cluster must be positive");
 
     const bool rotation_crystallography = opt.wedge_deg.has_value();
+    const bool selection = opt.selection;
 
     size_t nrefl = 0;
     for (const auto &i: observations)
@@ -744,7 +852,8 @@ ScaleMergeResult ScaleAndMergeReflectionsCeres(const std::vector<std::vector<Ref
         }
     }
 
-    scale(opt, g, mosaicity, R_sq, image_slot_used, rotation_crystallography, nhkl, obs);
+    std::cout << "Now scale the reflections: " << nrefl << std::endl;
+    scale(opt, g, mosaicity, R_sq, image_slot_used, rotation_crystallography, nhkl, obs, selection);
 
     ScaleMergeResult out;
 

image_analysis/scale_merge/ScaleAndMerge.h

@@ -50,6 +50,8 @@ struct ScaleMergeOptions {
 
     bool refine_wedge = false;
 
+    bool selection = true;
+
     enum class PartialityModel {Fixed, Rotation, Unity, Still} partiality_model = PartialityModel::Fixed;
 };
 

tools/jfjoch_process.cpp

@@ -49,8 +49,9 @@ void print_usage(Logger &logger) {
     logger.Info("   -A         Anomalous mode (don't merge Friedel pairs)");
     logger.Info("   -C<cell>   Fix reference unit cell: -C\"a,b,c,alpha,beta,gamma\" (comma-separated, no spaces; quotes optional)");
     logger.Info("   -c<num>    Max spot count (default: 250)");
-    logger.Info("   -W         HDF5 file with analysis results is written");
+    logger.Info("   -W<txt>    HDF5 file with analysis results is written. 'l' or 'light' deactivates image-output");
     logger.Info("   -T<num>    Noise sigma level for spot finding (default: 3.0)");
+    logger.Info("   -a         Use all reflection for scaling without filtering (default: false)");
 }
 
 void trim_in_place(std::string& t) {
@@ -163,6 +164,8 @@ int main(int argc, char **argv) {
     std::optional<int> space_group_number;
     std::optional<UnitCell> fixed_reference_unit_cell;
     bool write_output = false;
+    bool write_output_noimage = false;
+    bool filtering = true;
     std::optional<int64_t> max_spot_count_override;
     float sigma_spot_finding = 3.0;
     std::optional<float> merging_threshold;
@@ -179,7 +182,7 @@ int main(int argc, char **argv) {
     }
 
     int opt;
-    while ((opt = getopt(argc, argv, "o:N:s:e:vc:R::FX:xd:S:M::P:AD:C:T:W")) != -1) {
+    while ((opt = getopt(argc, argv, "o:N:s:e:vc:R::FX:xd:S:M::P:AD:C:T:W:a")) != -1) {
         switch (opt) {
             case 'o':
                 output_prefix = optarg;
@@ -195,6 +198,10 @@ int main(int argc, char **argv) {
                 break;
             case 'W':
                 write_output = true;
+                if (strcmp(optarg, "light") == 0 || strcmp(optarg, "l") == 0) {
+                    write_output_noimage = true;
+                    logger.Warning("Image data will not be saved.");
+                }
                 break;
             case 'v':
                 verbose = true;
@@ -237,6 +244,9 @@ int main(int argc, char **argv) {
             case 'x':
                 refine_beam_center = false;
                 break;
+            case 'a':
+                filtering = false;
+                break;
             case 'D':
                 d_min_scale_merge = atof(optarg);
                 logger.Info("High resolution limit for scaling/merging set to {:.2f} A", d_min_spot_finding);
@@ -334,8 +344,11 @@ int main(int argc, char **argv) {
     experiment.OverwriteExistingFiles(true);
     experiment.PolarizationFactor(0.99);
 
-    if (fixed_reference_unit_cell.has_value())
+    if (fixed_reference_unit_cell.has_value()) {
         experiment.SetUnitCell(*fixed_reference_unit_cell);
+    } else {
+        experiment.SetUnitCell({});
+    }
 
     if (max_spot_count_override.has_value()) {
         experiment.MaxSpotCount(max_spot_count_override.value());
@@ -346,6 +359,8 @@ int main(int argc, char **argv) {
     IndexingSettings indexing_settings;
     indexing_settings.Algorithm(indexing_algorithm);
     indexing_settings.RotationIndexing(rotation_indexing);
+    if (rotation_indexing)
+        logger.Info("Rotation indexing is activated.");
     if (rotation_indexing_range.has_value())
         indexing_settings.RotationIndexingMinAngularRange_deg(rotation_indexing_range.value());
 
@@ -355,6 +370,22 @@ int main(int argc, char **argv) {
         indexing_settings.GeomRefinementAlgorithm(GeomRefinementAlgorithmEnum::None);
     experiment.ImportIndexingSettings(indexing_settings);
 
+    switch (experiment.GetIndexingAlgorithm()) {
+        case IndexingAlgorithmEnum::FFBIDX:
+            logger.Info("Indexer used: FFBIDX");
+            break;
+        case IndexingAlgorithmEnum::FFTW:
+            logger.Info("Indexer used: FFTW");
+            break;
+        case IndexingAlgorithmEnum::FFT:
+            logger.Info("Indexer used: FFT (CUDA)");
+            break;
+        case IndexingAlgorithmEnum::None:
+            logger.Warning("Indexer not defined!");
+            return 0;
+        default: ;
+    }
+
     SpotFindingSettings spot_settings;
     spot_settings.enable = true;
     spot_settings.indexing = true;
@@ -435,6 +466,10 @@ int main(int argc, char **argv) {
         compressed_buffer.resize(MaxCompressedSize(experiment.GetCompressionAlgorithm(),
                                                    experiment.GetPixelsNum(),
                                                    experiment.GetByteDepthImage()));
+        auto size = compressor.Compress(compressed_buffer.data(),
+                            compressed_buffer.data(),
+                            experiment.GetPixelsNum(),
+                            sizeof(uint8_t));
 
         // Thread-local analysis resources
         MXAnalysisWithoutFPGA analysis(experiment, mapping, pixel_mask, indexer);
@@ -485,10 +520,11 @@ int main(int argc, char **argv) {
             auto image_end_time = std::chrono::high_resolution_clock::now();
             std::chrono::duration<float> image_duration = image_end_time - image_start_time;
 
-            auto size = compressor.Compress(compressed_buffer.data(),
-                                            img->Image().data(),
-                                            experiment.GetPixelsNum(),
-                                            sizeof(int32_t));
+            if (!write_output_noimage)
+                size = compressor.Compress(compressed_buffer.data(),
+                                                img->Image().data(),
+                                                experiment.GetPixelsNum(),
+                                                sizeof(int32_t));
 
             msg.image = CompressedImage(compressed_buffer.data(),
                                         size, experiment.GetXPixelsNum(),
@@ -521,7 +557,7 @@ int main(int argc, char **argv) {
             }
 
             // Progress log
-            if (current_idx_offset > 0 && current_idx_offset % 100 == 0) {
+            if ((current_idx_offset > 0 && (current_idx_offset+1) % 100 == 0) || image_idx == end_image - 1) {
                 std::optional<float> indexing_rate;
                 {
                     std::lock_guard<std::mutex> lock(plots_mutex);
@@ -530,11 +566,11 @@ int main(int argc, char **argv) {
 
                 if (indexing_rate.has_value()) {
                     logger.Info("Processed {} / {} images (indexing rate {:.1f}%)",
-                                current_idx_offset, images_to_process,
+                                current_idx_offset+1, images_to_process,
                                 indexing_rate.value() * 100.0f);
                 } else {
                     logger.Info("Processed {} / {} images (indexing rate N/A)",
-                                current_idx_offset, images_to_process);
+                                current_idx_offset+1, images_to_process);
                 }
                 if (image_idx == end_image - 1) {
                     if (!merging_threshold.has_value()) merging_threshold = 0.5f;
@@ -596,7 +632,6 @@ int main(int argc, char **argv) {
         logger.Info("Rotation Indexing found lattice");
     }
 
-    // --- Optional: run scaling (mosaicity refinement) on accumulated reflections ---
     // --- Optional: run scaling (mosaicity refinement) on accumulated reflections ---
     if (run_scaling) {
         logger.Info("Running scaling (mosaicity refinement) ...");
@@ -607,6 +642,7 @@ int main(int argc, char **argv) {
         scale_opts.max_solver_time_s = 240.0; // generous cutoff for now
         scale_opts.merge_friedel = !anomalous_mode;
         scale_opts.d_min_limit_A = d_min_scale_merge.value_or(0.0);
+        scale_opts.selection = filtering;
 
         const bool fixed_space_group = space_group || experiment.GetGemmiSpaceGroup().has_value();
 
@@ -728,6 +764,8 @@ int main(int argc, char **argv) {
                         cif_meta.unit_cell = rotation_indexer_ret->lattice.GetUnitCell();
                     } else if (experiment.GetUnitCell().has_value()) {
                         cif_meta.unit_cell = experiment.GetUnitCell().value();
+                    } else {
+                        logger.Warning("No UnitCell output");
                     }
 
                     if (scale_opts.space_group.has_value()) {