From bb40c9518bb1aec3e4e92928394cfc834696482a Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sun, 17 May 2026 14:44:33 +0200 Subject: [PATCH] Improve handling/calculating merge mask --- image_analysis/IndexAndRefine.cpp | 18 ++++++--------- image_analysis/IndexAndRefine.h | 3 ++- image_analysis/scale_merge/Merge.cpp | 17 +++++++++----- image_analysis/scale_merge/Merge.h | 2 +- tools/jfjoch_process.cpp | 33 ++++++++++++++++++---------- 5 files changed, 42 insertions(+), 31 deletions(-) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index fba19fb8..94abb0eb 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -304,15 +304,11 @@ ScalingResult IndexAndRefine::ScaleAllImages(const std::vector return result; } -MergeResult IndexAndRefine::Merge(bool calc_statistics, bool apply_cc_limit) const { - MergeResult out; - - std::vector merge_mask(reflections.size(), 1); - if (apply_cc_limit) - CalcMergeMask(experiment, scale_cc, merge_mask); - - out.merged = MergeAll(experiment, reflections, merge_mask); - if (calc_statistics) - out.statistics = MergeStats(experiment, out.merged, reflections, merge_mask); - return out; +const std::vector > &IndexAndRefine::GetReflections() const { + return reflections; } + +const std::vector &IndexAndRefine::GetImageCC() const { + return scale_cc; +} + diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index 6dc61b25..0e5553e1 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -68,7 +68,8 @@ public: IndexAndRefine& ReferenceIntensities(std::vector &reference); ScalingResult ScaleAllImages(const std::vector &reference, size_t nthreads = 0); - MergeResult Merge(bool statistics, bool apply_cc_limit) const; std::optional Finalize(); + const std::vector> &GetReflections() const; + const std::vector &GetImageCC() const; }; diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index bca14b14..9cf01704 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -16,19 +16,24 @@ // Need to use consistent lattice to calc resolution (which would also allow to calculate completeness) // But this is not possible for current still workflow (while fine for rotation workflow) -void CalcMergeMask(const DiffractionExperiment &x, +size_t CalcMergeMask(const DiffractionExperiment &x, const std::vector &scale_cc, std::vector &result_mask) { if (scale_cc.size() != result_mask.size()) throw JFJochException(JFJochExceptionCategory::InputParameterInvalid, "Scale CC size mismatch"); + size_t n_rejected = 0; auto min_cc = x.GetScalingSettings().GetMinCCForImage(); - if (min_cc <= 0.0) - return; - - for (int i = 0; i < result_mask.size(); ++i) - result_mask[i] = result_mask[i] && (scale_cc[i] < min_cc); + if (min_cc > 0.0) { + for (int i = 0; i < result_mask.size(); ++i) { + if (result_mask[i] && (scale_cc[i] < min_cc)) { + result_mask[i] = 0; + ++n_rejected; + } + } + } + return n_rejected; } std::vector MergeAll(const DiffractionExperiment &x, diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 3f2988fc..cd39b596 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -32,7 +32,7 @@ struct MergeResult { MergeStatistics statistics; }; -void CalcMergeMask(const DiffractionExperiment &x, +size_t CalcMergeMask(const DiffractionExperiment &x, const std::vector &scale_cc, std::vector &result_mask); diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index f029a06b..bf8cefab 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -604,6 +604,8 @@ int main(int argc, char **argv) { if (run_scaling || !reference_data.empty()) { logger.Info("Running scaling (mosaicity refinement) ..."); + std::vector scale_cc; + if (reference_data.empty()) { // If reference data are given, there is live scaling (no need to go again) ScalingResult scale_result(0); @@ -611,14 +613,15 @@ int main(int argc, char **argv) { auto scale_start = std::chrono::steady_clock::now(); for (int i = 0; i < scaling_iter; i++) { auto iter_start = std::chrono::steady_clock::now(); - auto merge_result = indexer.Merge(false, false); - scale_result = indexer.ScaleAllImages(merge_result.merged); + auto merge_result = MergeAll(experiment, indexer.GetReflections()); + scale_result = indexer.ScaleAllImages(merge_result); scale_result.SaveToFile(output_prefix + "_iter" + std::to_string(i) + "_scale.dat"); auto iter_end = std::chrono::steady_clock::now(); double iter_time = std::chrono::duration(iter_end - iter_start).count(); logger.Info("Scaling iteration {} took {:.3f} seconds", i, iter_time); } + scale_cc = scale_result.image_cc; end_msg.image_scale_factor = scale_result.image_scale_g; end_msg.image_scale_cc = scale_result.image_cc; end_msg.image_scale_mosaicity = scale_result.mosaicity_deg; @@ -627,15 +630,23 @@ int main(int argc, char **argv) { auto scale_end = std::chrono::steady_clock::now(); double scale_time = std::chrono::duration(scale_end - scale_start).count(); logger.Info("Scaling completed in {:.2f} s", scale_time); - } + } else + scale_cc = indexer.GetImageCC(); auto merge_start = std::chrono::steady_clock::now(); - auto merge_result = indexer.Merge(true, true); + + std::vector merge_mask(images_to_process, 1); + auto rejected = CalcMergeMask(experiment, scale_cc, merge_mask); + logger.Info("Rejected {} images due to low CC with reference", rejected); + + auto merged_reflections = MergeAll(experiment, indexer.GetReflections(), merge_mask); + auto merged_statistics = MergeStats(experiment, merged_reflections, indexer.GetReflections(), merge_mask); + auto merge_end = std::chrono::steady_clock::now(); double merge_time = std::chrono::duration(merge_end - merge_start).count(); logger.Info("Merge completed in {:.2f} s ({} unique reflections)", merge_time, - merge_result.merged.size()); + merged_reflections.size()); const bool fixed_space_group = space_group || experiment.GetGemmiSpaceGroup().has_value(); @@ -653,7 +664,7 @@ int main(int argc, char **argv) { sg_opts.min_total_compared = 100; sg_opts.test_systematic_absences = true; - const auto sg_search = SearchSpaceGroup(merge_result.merged, sg_opts); + const auto sg_search = SearchSpaceGroup(merged_reflections, sg_opts); logger.Info(""); { @@ -686,9 +697,7 @@ int main(int argc, char **argv) { } */ } - // Print resolution-shell statistics table - const auto &stats = merge_result.statistics; - stats.Print(logger); + merged_statistics.Print(logger); { const std::string hkl_path = output_prefix + "_intensities.hkl"; @@ -696,13 +705,13 @@ int main(int argc, char **argv) { if (!hkl_file) { logger.Error("Cannot open {} for writing", hkl_path); } else { - for (const auto &r: merge_result.merged) { + for (const auto &r: merged_reflections) { hkl_file << r.h << " " << r.k << " " << r.l << " " << r.I << " " << r.sigma << "\n"; } hkl_file.close(); - logger.Info("Wrote {} reflections to {}", merge_result.merged.size(), hkl_path); + logger.Info("Wrote {} reflections to {}", merged_reflections.size(), hkl_path); } } @@ -717,7 +726,7 @@ int main(int argc, char **argv) { cif_meta.unit_cell = rotation_indexer_ret->lattice.GetUnitCell(); cif_meta.data_block_name = output_prefix; - WriteMmcifReflections(cif_path, merge_result.merged, cif_meta); + WriteMmcifReflections(cif_path, merged_reflections, cif_meta); logger.Info("Wrote mmCIF reflections to {}", cif_path); } catch (const std::exception &e) { logger.Error("Failed to write mmCIF: {}", e.what());