From 1bf70dcae732ec9d3b69d513fe4318a8d5ac1b2b Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Sun, 17 May 2026 14:32:39 +0200 Subject: [PATCH] Merge: Some clean-up of CC merge mask + add unit cell handling in IndexAndRefine --- image_analysis/IndexAndRefine.cpp | 23 +++++++++++++---------- image_analysis/IndexAndRefine.h | 3 ++- image_analysis/scale_merge/Merge.cpp | 16 ++++++++++------ image_analysis/scale_merge/Merge.h | 5 +++-- tools/jfjoch_scale.cpp | 4 +++- 5 files changed, 31 insertions(+), 20 deletions(-) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index 275e28d3..fba19fb8 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -20,7 +20,8 @@ IndexAndRefine::IndexAndRefine(const DiffractionExperiment &x, IndexerThreadPool rotation_indexer = std::make_unique(x, *indexer); reflections.resize(x.GetImageNum()); mosaicity.resize(x.GetImageNum(), NAN); - merge_mask.resize(x.GetImageNum(), 0); + scale_cc.resize(x.GetImageNum(), 0); + unit_cells.resize(x.GetImageNum()); } IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetryRotation(DataMessage &msg) { @@ -256,12 +257,13 @@ void IndexAndRefine::ProcessImage(DataMessage &msg, if (experiment.GetIndexingSettings().GetGeomRefinementAlgorithm() != GeomRefinementAlgorithmEnum::None) RefineGeometryIfNeeded(msg, outcome); - if (!outcome.lattice_candidate) + if (!outcome.lattice_candidate.has_value()) return; if (!AnalyzeIndexing(msg, outcome.experiment, *outcome.lattice_candidate)) return; + unit_cells[msg.number] = outcome.lattice_candidate->GetUnitCell(); msg.lattice_type = outcome.symmetry; if (spot_finding_settings.quick_integration) @@ -289,10 +291,7 @@ void IndexAndRefine::ScaleImage(DataMessage &msg) { msg.image_scale_factor = scaling_result.G; msg.image_scale_mosaicity = scaling_result.mos; msg.image_scale_cc = scaling_result.cc; - if (scaling_result.cc >= experiment.GetScalingSettings().GetMinCCForImage()) - merge_mask[msg.number] = 1; - else - merge_mask[msg.number] = 0; + scale_cc[msg.number] = scaling_result.cc; auto scaling_end_time = std::chrono::steady_clock::now(); msg.image_scale_time_s = std::chrono::duration(scaling_end_time - scaling_start_time).count(); @@ -301,15 +300,19 @@ void IndexAndRefine::ScaleImage(DataMessage &msg) { ScalingResult IndexAndRefine::ScaleAllImages(const std::vector &reference, size_t nthreads) { ScaleOnTheFly scaling(experiment, reference); auto result = scaling.Scale(reflections, mosaicity, nthreads); - merge_mask = CalcMergeMask(experiment, result); + scale_cc = result.image_cc; return result; } MergeResult IndexAndRefine::Merge(bool calc_statistics, bool apply_cc_limit) const { MergeResult out; - const auto& mask = apply_cc_limit ? merge_mask : std::vector{}; - out.merged = MergeAll(experiment, reflections, mask); + + 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, mask); + out.statistics = MergeStats(experiment, out.merged, reflections, merge_mask); return out; } diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index ad1734f5..6dc61b25 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -48,7 +48,8 @@ class IndexAndRefine { mutable std::mutex reflections_mutex; std::vector> reflections; std::vector mosaicity; - std::vector merge_mask; + std::vector scale_cc; + std::vector > unit_cells; IndexingOutcome DetermineLatticeAndSymmetryRotation(DataMessage &msg); IndexingOutcome DetermineLatticeAndSymmetry(DataMessage &msg); diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index f7c55350..bca14b14 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -16,15 +16,19 @@ // 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) -std::vector CalcMergeMask(const DiffractionExperiment &x, - const ScalingResult &result) { - std::vector ret(result.image_cc.size(), 0); +void 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"); auto min_cc = x.GetScalingSettings().GetMinCCForImage(); - for (int i = 0 ; i < result.image_cc.size(); ++i) - ret[i] = (result.image_cc[i] >= min_cc) ? 1 : 0; + if (min_cc <= 0.0) + return; - return ret; + for (int i = 0; i < result_mask.size(); ++i) + result_mask[i] = result_mask[i] && (scale_cc[i] < min_cc); } std::vector MergeAll(const DiffractionExperiment &x, diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index e08bbeb6..3f2988fc 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -32,8 +32,9 @@ struct MergeResult { MergeStatistics statistics; }; -std::vector CalcMergeMask(const DiffractionExperiment &x, - const ScalingResult &result); +void CalcMergeMask(const DiffractionExperiment &x, + const std::vector &scale_cc, + std::vector &result_mask); std::vector MergeAll(const DiffractionExperiment &x, const std::vector > &reflections, diff --git a/tools/jfjoch_scale.cpp b/tools/jfjoch_scale.cpp index 42c726f9..e5db9be3 100644 --- a/tools/jfjoch_scale.cpp +++ b/tools/jfjoch_scale.cpp @@ -271,9 +271,11 @@ int main(int argc, char **argv) { double scale_time = std::chrono::duration(scale_end - scale_start).count(); logger.Info("Scaling completed in {:.2f} s", scale_time); + std::vector merge_mask(images_to_process, 1); + auto merge_start = std::chrono::steady_clock::now(); - auto merge_mask = CalcMergeMask(experiment, scale_result); + CalcMergeMask(experiment, scale_result.image_cc, merge_mask); auto merge_result = MergeAll(experiment, reflections, merge_mask); auto merge_stats = MergeStats(experiment, merge_result, reflections, merge_mask);