From 50eeba40059b37a22d8e2cec4138a379f9c98eaa Mon Sep 17 00:00:00 2001 From: Mazzoleni Alice Francesca Date: Wed, 2 Apr 2025 17:58:26 +0200 Subject: [PATCH] restructured GainMap to have own class and generalized --- CMakeLists.txt | 1 + include/aare/ClusterFile.hpp | 30 ++++++++++++----- include/aare/ClusterVector.hpp | 24 +------------- include/aare/GainMap.hpp | 59 ++++++++++++++++++++++++++++++++++ src/ClusterVector.test.cpp | 46 ++++++++++++++++++++++++++ 5 files changed, 129 insertions(+), 31 deletions(-) create mode 100644 include/aare/GainMap.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 0ab1e73..b02303c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -344,6 +344,7 @@ set(PUBLICHEADERS include/aare/Fit.hpp include/aare/FileInterface.hpp include/aare/Frame.hpp + include/aare/GainMap.hpp include/aare/geo_helpers.hpp include/aare/NDArray.hpp include/aare/NDView.hpp diff --git a/include/aare/ClusterFile.hpp b/include/aare/ClusterFile.hpp index 9c43326..eb6cc86 100644 --- a/include/aare/ClusterFile.hpp +++ b/include/aare/ClusterFile.hpp @@ -2,6 +2,7 @@ #include "aare/Cluster.hpp" #include "aare/ClusterVector.hpp" +#include "aare/GainMap.hpp" #include "aare/NDArray.hpp" #include "aare/defs.hpp" #include @@ -44,9 +45,8 @@ class ClusterFile { std::optional m_roi; /*Region of interest, will be applied if set*/ std::optional> m_noise_map; /*Noise map to cut photons, will be applied if set*/ - std::optional> - m_gain_map; /*Gain map to apply to the clusters, will be applied if - set*/ + std::optional m_gain_map; /*Gain map to apply to the clusters, will + be applied if set*/ public: /** @@ -107,6 +107,10 @@ class ClusterFile { */ void set_gain_map(const NDView gain_map); + void set_gain_map(const GainMap &gain_map); + + void set_gain_map(const GainMap &&gain_map); + /** * @brief Close the file. If not closed the file will be closed in the * destructor @@ -175,7 +179,17 @@ void ClusterFile::set_noise_map( template void ClusterFile::set_gain_map( const NDView gain_map) { - m_gain_map = NDArray(gain_map); + m_gain_map = GainMap(gain_map); +} + +template +void ClusterFile::set_gain_map(const GainMap &gain_map) { + m_gain_map = gain_map; +} + +template +void ClusterFile::set_gain_map(const GainMap &&gain_map) { + m_gain_map = gain_map; } // TODO generally supported for all clsuter types @@ -263,7 +277,7 @@ ClusterFile::read_clusters_without_cut(size_t n_clusters) { // No new allocation, only change bounds. clusters.resize(nph_read); if (m_gain_map) - clusters.apply_gain_map(m_gain_map->view()); + m_gain_map->apply_gain_map(clusters); return clusters; } @@ -312,7 +326,7 @@ ClusterFile::read_clusters_with_cut(size_t n_clusters) { } } if (m_gain_map) - clusters.apply_gain_map(m_gain_map->view()); + m_gain_map->apply_gain_map(clusters); return clusters; } @@ -370,7 +384,7 @@ ClusterFile::read_frame_without_cut() { } clusters.resize(n_clusters); if (m_gain_map) - clusters.apply_gain_map(m_gain_map->view()); + m_gain_map->apply_gain_map(clusters); return clusters; } @@ -403,7 +417,7 @@ ClusterFile::read_frame_with_cut() { } } if (m_gain_map) - clusters.apply_gain_map(m_gain_map->view()); + m_gain_map->apply_gain_map(clusters); return clusters; } diff --git a/include/aare/ClusterVector.hpp b/include/aare/ClusterVector.hpp index 30be5eb..ca2fd4d 100644 --- a/include/aare/ClusterVector.hpp +++ b/include/aare/ClusterVector.hpp @@ -32,7 +32,6 @@ class ClusterVector; // Forward declaration template class ClusterVector> { - using value_type = T; std::byte *m_data{}; size_t m_size{0}; @@ -49,6 +48,7 @@ class ClusterVector> { constexpr static char m_fmt_base[] = "=h:x:\nh:y:\n({},{}){}:data:"; public: + using value_type = T; using ClusterType = Cluster; /** @@ -237,28 +237,6 @@ class ClusterVector> { m_size = new_size; } - // TODO: Generalize !!!! Maybe move somewhere else - void apply_gain_map(const NDView gain_map) { - // in principle we need to know the size of the image for this lookup - // TODO! check orientations - std::array xcorr = {-1, 0, 1, -1, 0, 1, -1, 0, 1}; - std::array ycorr = {-1, -1, -1, 0, 0, 0, 1, 1, 1}; - for (size_t i = 0; i < m_size; i++) { - auto &cl = at(i); - - if (cl.x > 0 && cl.y > 0 && cl.x < gain_map.shape(1) - 1 && - cl.y < gain_map.shape(0) - 1) { - for (size_t j = 0; j < 9; j++) { - size_t x = cl.x + xcorr[j]; - size_t y = cl.y + ycorr[j]; - cl.data[j] = static_cast(cl.data[j] * gain_map(y, x)); - } - } else { - memset(cl.data, 0, 9 * sizeof(T)); // clear edge clusters - } - } - } - private: void allocate_buffer(size_t new_capacity) { size_t num_bytes = item_size() * new_capacity; diff --git a/include/aare/GainMap.hpp b/include/aare/GainMap.hpp new file mode 100644 index 0000000..9eb7c11 --- /dev/null +++ b/include/aare/GainMap.hpp @@ -0,0 +1,59 @@ +/************************************************ + * @file ApplyGainMap.hpp + * @short function to apply gain map of image size to a vector of clusters + ***********************************************/ + +#pragma once +#include "aare/Cluster.hpp" +#include "aare/ClusterVector.hpp" +#include "aare/NDArray.hpp" +#include "aare/NDView.hpp" +#include + +namespace aare { + +class GainMap { + + public: + explicit GainMap(const NDArray &gain_map) + : m_gain_map(gain_map) {}; + + explicit GainMap(const NDView gain_map) { + m_gain_map = NDArray(gain_map); + } + + template >> + void apply_gain_map(ClusterVector &clustervec) { + // in principle we need to know the size of the image for this lookup + // TODO! check orientations + size_t ClusterSizeX = clustervec.cluster_size_x(); + size_t ClusterSizeY = clustervec.cluster_size_y(); + + using T = typename ClusterVector::value_type; + + int64_t index_cluster_center_x = ClusterSizeX / 2; + int64_t index_cluster_center_y = ClusterSizeY / 2; + for (size_t i = 0; i < clustervec.size(); i++) { + auto &cl = clustervec.at(i); + + if (cl.x > 0 && cl.y > 0 && cl.x < m_gain_map.shape(1) - 1 && + cl.y < m_gain_map.shape(0) - 1) { + for (size_t j = 0; j < ClusterSizeX * ClusterSizeY; j++) { + size_t x = cl.x + j % ClusterSizeX - index_cluster_center_x; + size_t y = cl.y + j / ClusterSizeX - index_cluster_center_y; + cl.data[j] = static_cast(cl.data[j] * m_gain_map(y, x)); + } + } else { + memset(cl.data, 0, + ClusterSizeX * ClusterSizeY * + sizeof(T)); // clear edge clusters + } + } + } + + private: + NDArray m_gain_map{}; +}; + +} // end of namespace aare \ No newline at end of file diff --git a/src/ClusterVector.test.cpp b/src/ClusterVector.test.cpp index b58e88a..c354891 100644 --- a/src/ClusterVector.test.cpp +++ b/src/ClusterVector.test.cpp @@ -1,6 +1,7 @@ #include "aare/ClusterVector.hpp" #include +#include #include #include @@ -183,4 +184,49 @@ TEST_CASE("Concatenate two cluster vectors where we need to allocate", REQUIRE(ptr[2].y == 12); REQUIRE(ptr[3].x == 16); REQUIRE(ptr[3].y == 17); +} + +struct ClusterTestData { + int8_t ClusterSizeX; + int8_t ClusterSizeY; + std::vector index_map_x; + std::vector index_map_y; +}; + +TEST_CASE("Gain Map Calculation Index Map", "[.ClusterVector][.gain_map]") { + + auto clustertestdata = GENERATE( + ClusterTestData{3, + 3, + {-1, 0, 1, -1, 0, 1, -1, 0, 1}, + {-1, -1, -1, 0, 0, 0, 1, 1, 1}}, + ClusterTestData{ + 4, + 4, + {-2, -1, 0, 1, -2, -1, 0, 1, -2, -1, 0, 1, -2, -1, 0, 1}, + {-2, -2, -2, -2, -1, -1, -1, -1, 0, 0, 0, 0, 1, 1, 1, 1}}, + ClusterTestData{2, 2, {-1, 0, -1, 0}, {-1, -1, 0, 0}}, + ClusterTestData{5, + 5, + {-2, -1, 0, 1, 2, -2, -1, 0, 1, 2, -2, -1, 0, + 1, 2, -2, -1, 0, 1, 2, -2, -1, 0, 1, 2}, + {-2, -2, -2, -2, -2, -1, -1, -1, -1, -1, 0, 0, 0, + 0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2}}); + + int8_t ClusterSizeX = clustertestdata.ClusterSizeX; + int8_t ClusterSizeY = clustertestdata.ClusterSizeY; + + std::vector index_map_x(ClusterSizeX * ClusterSizeY); + std::vector index_map_y(ClusterSizeX * ClusterSizeY); + + int64_t index_cluster_center_x = ClusterSizeX / 2; + int64_t index_cluster_center_y = ClusterSizeY / 2; + + for (size_t j = 0; j < ClusterSizeX * ClusterSizeY; j++) { + index_map_x[j] = j % ClusterSizeX - index_cluster_center_x; + index_map_y[j] = j / ClusterSizeX - index_cluster_center_y; + } + + CHECK(index_map_x == clustertestdata.index_map_x); + CHECK(index_map_y == clustertestdata.index_map_y); } \ No newline at end of file