From 2397a108dd7185c8347667b26b9cbdfc3a4903d7 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 13 May 2026 15:23:32 +0200 Subject: [PATCH] ScaleOnTheFly: Add CC image/ref --- image_analysis/scale_merge/ScaleOnTheFly.cpp | 59 ++++++++++++++++++++ image_analysis/scale_merge/ScaleOnTheFly.h | 3 + image_analysis/scale_merge/ScalingResult.cpp | 11 ++-- image_analysis/scale_merge/ScalingResult.h | 3 +- tools/jfjoch_process.cpp | 16 +++--- 5 files changed, 78 insertions(+), 14 deletions(-) diff --git a/image_analysis/scale_merge/ScaleOnTheFly.cpp b/image_analysis/scale_merge/ScaleOnTheFly.cpp index d2057fbe0..b1e63ac4e 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.cpp +++ b/image_analysis/scale_merge/ScaleOnTheFly.cpp @@ -115,6 +115,57 @@ bool ScaleOnTheFly::Accept(const Reflection &r) { return true; } +std::pair ScaleOnTheFly::CalculateGlobalCC(const std::vector &reflections) const { + double sum_x = 0.0; + double sum_y = 0.0; + double sum_x2 = 0.0; + double sum_y2 = 0.0; + double sum_xy = 0.0; + size_t n = 0; + + for (const auto &r: reflections) { + if (!AcceptReflection(r, s.GetHighResolutionLimit_A())) + continue; + if (r.partiality < s.GetMinPartiality()) + continue; + if (!std::isfinite(r.I) || !std::isfinite(r.scaling_correction) || r.scaling_correction <= 0.0f) + continue; + if (!std::isfinite(r.sigma) || r.sigma <= 0.0f) + continue; + + const HKLKey key = hkl_key_generator(r); + const auto it = reference_data.find(key); + if (it == reference_data.end()) + continue; + + const double image_i = static_cast(r.I) * static_cast(r.scaling_correction); + const double ref_i = it->second; + + if (!std::isfinite(image_i) || !std::isfinite(ref_i)) + continue; + + sum_x += image_i; + sum_y += ref_i; + sum_x2 += image_i * image_i; + sum_y2 += ref_i * ref_i; + sum_xy += image_i * ref_i; + ++n; + } + + if (n < MIN_REFLECTIONS) + return {NAN, n}; + + const double nd = static_cast(n); + const double cov = sum_xy - sum_x * sum_y / nd; + const double var_x = sum_x2 - sum_x * sum_x / nd; + const double var_y = sum_y2 - sum_y * sum_y / nd; + + if (!(var_x > 0.0 && var_y > 0.0)) + return {NAN, n}; + + return {cov / std::sqrt(var_x * var_y), n}; +} + ScaleOnTheFlyResult ScaleOnTheFly::Scale(std::vector &reflections, std::optional mosaicity_deg) { auto start = std::chrono::steady_clock::now(); @@ -235,6 +286,10 @@ ScaleOnTheFlyResult ScaleOnTheFly::Scale(std::vector &reflections, s } } + const auto [cc, cc_n] = CalculateGlobalCC(reflections); + result.cc = cc; + result.cc_n = cc_n; + auto end = std::chrono::steady_clock::now(); result.time_s = std::chrono::duration(end - start).count(); return result; @@ -262,6 +317,8 @@ ScalingResult ScaleOnTheFly::Scale(std::vector > &reflec result.image_bfactor_Ang2[i] = local_result.B; result.image_scale_g[i] = local_result.G; result.rotation_wedge_deg[i] = local_result.wedge; + result.image_cc[i] = local_result.cc; + result.image_cc_n[i] = local_result.cc_n; } } else { auto local_nthreads = std::min(nthreads, reflections.size()); @@ -285,6 +342,8 @@ ScalingResult ScaleOnTheFly::Scale(std::vector > &reflec result.image_bfactor_Ang2[i] = local_result.B; result.image_scale_g[i] = local_result.G; result.rotation_wedge_deg[i] = local_result.wedge; + result.image_cc[i] = local_result.cc; + result.image_cc_n[i] = local_result.cc_n; i = curr_image.fetch_add(1); } })); diff --git a/image_analysis/scale_merge/ScaleOnTheFly.h b/image_analysis/scale_merge/ScaleOnTheFly.h index c491c4a51..fac68f3bf 100644 --- a/image_analysis/scale_merge/ScaleOnTheFly.h +++ b/image_analysis/scale_merge/ScaleOnTheFly.h @@ -16,6 +16,8 @@ struct ScaleOnTheFlyResult { double G = 1.0; double mos = 0.1; double wedge = 0.1; + double cc = NAN; + size_t cc_n = 0; float time_s = 0.0; bool succesful = false; }; @@ -32,6 +34,7 @@ class ScaleOnTheFly { std::map reference_data; bool Accept(const Reflection &r); + [[nodiscard]] std::pair CalculateGlobalCC(const std::vector &reflections) const; public: ScaleOnTheFly(const std::vector &ref, const DiffractionExperiment &x); ScaleOnTheFlyResult Scale(std::vector &r, std::optional mosaicity_deg); diff --git a/image_analysis/scale_merge/ScalingResult.cpp b/image_analysis/scale_merge/ScalingResult.cpp index 8c683da99..a5dc4aaea 100644 --- a/image_analysis/scale_merge/ScalingResult.cpp +++ b/image_analysis/scale_merge/ScalingResult.cpp @@ -11,25 +11,26 @@ ScalingResult::ScalingResult(size_t n) : image_scale_g(n, NAN), mosaicity_deg(n, NAN), image_bfactor_Ang2(n, NAN), - rotation_wedge_deg(n, NAN) { -} + rotation_wedge_deg(n, NAN), + image_cc(n, NAN), + image_cc_n(n, 0) {} void ScalingResult::SaveToFile(const std::string &filename) { const std::string img_path = filename + "_image.dat"; - std::ofstream img_file(img_path); + std::ofstream img_file(img_path, std::ofstream::out | std::ofstream::trunc); if (!img_file) { throw JFJochException(JFJochExceptionCategory::FileWriteError , "Cannot open {} for writing"); } - img_file << "# image_id G B mosaicity_deg wedge_deg\n"; - for (size_t i = 0; i < image_scale_g.size(); ++i) { img_file << i << " " << image_scale_g[i] << " " << image_bfactor_Ang2[i] << " " << mosaicity_deg[i] << " " << rotation_wedge_deg[i] + << " " << image_cc[i] + << " " << image_cc_n[i] << "\n"; } diff --git a/image_analysis/scale_merge/ScalingResult.h b/image_analysis/scale_merge/ScalingResult.h index 880f5c55b..91bc8c2a1 100644 --- a/image_analysis/scale_merge/ScalingResult.h +++ b/image_analysis/scale_merge/ScalingResult.h @@ -11,7 +11,8 @@ struct ScalingResult { std::vector mosaicity_deg; std::vector image_bfactor_Ang2; std::vector rotation_wedge_deg; - + std::vector image_cc; + std::vector image_cc_n; explicit ScalingResult(size_t n); void SaveToFile(const std::string &filename); }; diff --git a/tools/jfjoch_process.cpp b/tools/jfjoch_process.cpp index 9c2d57690..5eceeed09 100644 --- a/tools/jfjoch_process.cpp +++ b/tools/jfjoch_process.cpp @@ -587,20 +587,20 @@ int main(int argc, char **argv) { const bool fixed_space_group = space_group || experiment.GetGemmiSpaceGroup().has_value(); + ScalingResult scale_result(0); + auto scale_start = std::chrono::steady_clock::now(); for (int i = 0; i < 3; i++) { auto iter_start = std::chrono::steady_clock::now(); if (reference_data.empty()) { auto merge_result = indexer.Merge(false); - auto scale_result = indexer.ScaleAllImages(merge_result.merged); - end_msg.image_scale_factor = scale_result.image_scale_g; - scale_result.SaveToFile(output_prefix + "_iter" + std::to_string(i) + "_scale.dat"); - } else { - auto scale_result = indexer.ScaleAllImages(reference_data); - end_msg.image_scale_factor = scale_result.image_scale_g; - scale_result.SaveToFile(output_prefix + "_iter" + std::to_string(i) + "_scale.dat"); - } + scale_result = indexer.ScaleAllImages(merge_result.merged); + } else + scale_result = indexer.ScaleAllImages(reference_data); + + end_msg.image_scale_factor = scale_result.image_scale_g; + 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();