ScaleOnTheFly: Add CC image/ref

image_analysis/scale_merge/ScaleOnTheFly.cpp

@@ -115,6 +115,57 @@ bool ScaleOnTheFly::Accept(const Reflection &r) {
     return true;
 }
 
+std::pair<double, size_t> ScaleOnTheFly::CalculateGlobalCC(const std::vector<Reflection> &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<double>(r.I) * static_cast<double>(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<double>(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<Reflection> &reflections, std::optional<double> mosaicity_deg) {
     auto start = std::chrono::steady_clock::now();
 
@@ -235,6 +286,10 @@ ScaleOnTheFlyResult ScaleOnTheFly::Scale(std::vector<Reflection> &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<float>(end - start).count();
     return result;
@@ -262,6 +317,8 @@ ScalingResult ScaleOnTheFly::Scale(std::vector<std::vector<Reflection> > &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<std::vector<Reflection> > &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);
                 }
             }));

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<HKLKey, double> reference_data;
 
     bool Accept(const Reflection &r);
+    [[nodiscard]] std::pair<double, size_t> CalculateGlobalCC(const std::vector<Reflection> &reflections) const;
 public:
     ScaleOnTheFly(const std::vector<MergedReflection> &ref, const DiffractionExperiment &x);
     ScaleOnTheFlyResult Scale(std::vector<Reflection> &r, std::optional<double> mosaicity_deg);

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";
     }
 

image_analysis/scale_merge/ScalingResult.h

@@ -11,7 +11,8 @@ struct ScalingResult {
     std::vector<float> mosaicity_deg;
     std::vector<float> image_bfactor_Ang2;
     std::vector<float> rotation_wedge_deg;
-
+    std::vector<float> image_cc;
+    std::vector<int> image_cc_n;
     explicit ScalingResult(size_t n);
     void SaveToFile(const std::string &filename);
 };

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<double>(iter_end - iter_start).count();