From c90e4facce452cee96fe7024dc4317c580cf4b79 Mon Sep 17 00:00:00 2001 From: leonarski_f Date: Mon, 18 May 2026 20:23:13 +0200 Subject: [PATCH] CorrelationCoefficient: Separate class --- common/CMakeLists.txt | 2 ++ common/CorrelationCoefficient.cpp | 33 ++++++++++++++++++++ common/CorrelationCoefficient.h | 17 +++++++++++ image_analysis/scale_merge/Merge.cpp | 45 ++++------------------------ image_analysis/scale_merge/Merge.h | 1 - tests/CCTest.cpp | 22 ++++++++++++++ tests/CMakeLists.txt | 1 + 7 files changed, 81 insertions(+), 40 deletions(-) create mode 100644 common/CorrelationCoefficient.cpp create mode 100644 common/CorrelationCoefficient.h create mode 100644 tests/CCTest.cpp diff --git a/common/CMakeLists.txt b/common/CMakeLists.txt index 0f5b5574..14c2d364 100644 --- a/common/CMakeLists.txt +++ b/common/CMakeLists.txt @@ -130,6 +130,8 @@ ADD_LIBRARY(JFJochCommon STATIC ScalingSettings.h UnitCell.cpp UnitCell.h + CorrelationCoefficient.cpp + CorrelationCoefficient.h ) TARGET_LINK_LIBRARIES(JFJochCommon JFJochLogger Compression JFCalibration gemmi Threads::Threads -lrt ) diff --git a/common/CorrelationCoefficient.cpp b/common/CorrelationCoefficient.cpp new file mode 100644 index 00000000..24925571 --- /dev/null +++ b/common/CorrelationCoefficient.cpp @@ -0,0 +1,33 @@ +// +// Created by leonarski_f on 18.05.2026. +// + +#include + +#include "CorrelationCoefficient.h" + +void CorrelationCoefficient::Add(double x, double y) { + sum_x += x; + sum_y += y; + sum_x2 += x * x; + sum_y2 += y * y; + sum_xy += x * y; + n_cc++; +} + +double CorrelationCoefficient::GetCC() const { + if (n_cc < 2) + return NAN; + const double n = static_cast(n_cc); + const double mean_x = sum_x / n; + const double mean_y = sum_y / n; + + const double cov = sum_xy / n - mean_x * mean_y; + const double var_x = sum_x2 / n - mean_x * mean_x; + const double var_y = sum_y2 / n - mean_y * mean_y; + + if (!(var_x > 0.0 && var_y > 0.0)) + return NAN; + + return cov / std::sqrt(var_x * var_y); +} diff --git a/common/CorrelationCoefficient.h b/common/CorrelationCoefficient.h new file mode 100644 index 00000000..71128422 --- /dev/null +++ b/common/CorrelationCoefficient.h @@ -0,0 +1,17 @@ +// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +class CorrelationCoefficient { + 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; + int n_cc = 0; +public: + void Add(double x, double y); + [[nodiscard]] double GetCC() const; +}; + diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index 6f85e66d..4bc2bec4 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -11,6 +11,7 @@ #include +#include "../../common/CorrelationCoefficient.h" #include "../../common/ResolutionShells.h" #include "HKLKey.h" @@ -170,39 +171,6 @@ std::vector MergeAll(const DiffractionExperiment &x, return out; } -struct CCAccum { - 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; - int n_cc = 0; - void Add(double x, double y) { - sum_x += x; - sum_y += y; - sum_x2 += x * x; - sum_y2 += y * y; - sum_xy += x * y; - n_cc++; - } - - double GetCC() const { - if (n_cc < 2) - return NAN; - const double n = static_cast(n_cc); - const double mean_x = sum_x / n; - const double mean_y = sum_y / n; - - const double cov = sum_xy / n - mean_x * mean_y; - const double var_x = sum_x2 / n - mean_x * mean_x; - const double var_y = sum_y2 / n - mean_y * mean_y; - - if (!(var_x > 0.0 && var_y > 0.0)) - return NAN; - - return cov / std::sqrt(var_x * var_y); - } -}; struct ShellAccum { int total_obs = 0; @@ -212,8 +180,8 @@ struct ShellAccum { double sum_i_over_sigma = 0.0; int n_i_over_sigma = 0; - CCAccum cc_half; - CCAccum cc_ref; + CorrelationCoefficient cc_half; + CorrelationCoefficient cc_ref; }; void CalcPossibleReflections(const DiffractionExperiment &x, @@ -303,8 +271,8 @@ MergeStatistics MergeStats(const DiffractionExperiment &x, CalcPossibleReflections(x, cell, d_min_pad, d_max_pad, shells, acc); - CCAccum cc_half_overall; - CCAccum cc_ref_overall; + CorrelationCoefficient cc_half_overall; + CorrelationCoefficient cc_ref_overall; for (const auto &m: merged) { const auto shell = shells.GetShell(m.d); @@ -326,6 +294,7 @@ MergeStatistics MergeStats(const DiffractionExperiment &x, cc_ref_overall.Add(m.I, ref_it->second); } } + if (std::isfinite(m.I_half[0]) && std::isfinite(m.I_half[1])) { acc[s].cc_half.Add(m.I_half[0], m.I_half[1]); cc_half_overall.Add(m.I_half[0], m.I_half[1]); @@ -342,10 +311,8 @@ MergeStatistics MergeStats(const DiffractionExperiment &x, for (const auto &r: reflections[i]) { if (key_generator.IsSystematicallyAbsent(r)) continue; - if (r.image_scale_corr <= 0.0 || !std::isfinite(r.image_scale_corr)) continue; - if (!AcceptReflection(r, d_min_limit_A)) continue; if (r.partiality < min_partiality) diff --git a/image_analysis/scale_merge/Merge.h b/image_analysis/scale_merge/Merge.h index 5cf6b000..529d8dba 100644 --- a/image_analysis/scale_merge/Merge.h +++ b/image_analysis/scale_merge/Merge.h @@ -10,7 +10,6 @@ #include "../../common/Reflection.h" struct MergeStatisticsShell { - float d_min = 0.0f; float d_max = 0.0f; float mean_one_over_d2 = 0; diff --git a/tests/CCTest.cpp b/tests/CCTest.cpp new file mode 100644 index 00000000..328f5147 --- /dev/null +++ b/tests/CCTest.cpp @@ -0,0 +1,22 @@ +// SPDX-FileCopyrightText: 2024 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include +#include "../common/CorrelationCoefficient.h" + +TEST_CASE("CorrelationCoefficient") { + CorrelationCoefficient cc; + CHECK(std::isnan(cc.GetCC())); + + cc.Add(100.0, 500.0); + CHECK(std::isnan(cc.GetCC())); + + cc.Add(200.0, 1000.0); + CHECK(cc.GetCC() == Catch::Approx(1.0)); + + cc.Add(300.0, 1500.0); + CHECK(cc.GetCC() == Catch::Approx(1.0)); + + cc.Add(400.0, 2000.0); + CHECK(cc.GetCC() == Catch::Approx(1.0)); +} diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index b38c8e97..951350cd 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -69,6 +69,7 @@ ADD_EXECUTABLE(jfjoch_test XDSPluginTest.cpp MergeScaleTest.cpp UnitCellTest.cpp + CCTest.cpp ) target_link_libraries(jfjoch_test Catch2WithMain JFJochBroker JFJochReceiver JFJochReader JFJochWriter