restructured GainMap to have own class and generalized

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

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 <filesystem>
@@ -44,9 +45,8 @@ class ClusterFile {
     std::optional<ROI> m_roi; /*Region of interest, will be applied if set*/
     std::optional<NDArray<int32_t, 2>>
         m_noise_map; /*Noise map to cut photons, will be applied if set*/
-    std::optional<NDArray<double, 2>>
-        m_gain_map; /*Gain map to apply to the clusters, will be applied if
-                       set*/
+    std::optional<GainMap> 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<double, 2> 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<ClusterType, Enable>::set_noise_map(
 template <typename ClusterType, typename Enable>
 void ClusterFile<ClusterType, Enable>::set_gain_map(
     const NDView<double, 2> gain_map) {
-    m_gain_map = NDArray<double, 2>(gain_map);
+    m_gain_map = GainMap(gain_map);
+}
+
+template <typename ClusterType, typename Enable>
+void ClusterFile<ClusterType, Enable>::set_gain_map(const GainMap &gain_map) {
+    m_gain_map = gain_map;
+}
+
+template <typename ClusterType, typename Enable>
+void ClusterFile<ClusterType, Enable>::set_gain_map(const GainMap &&gain_map) {
+    m_gain_map = gain_map;
 }
 
 // TODO generally supported for all clsuter types
@@ -263,7 +277,7 @@ ClusterFile<ClusterType, Enable>::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<ClusterType, Enable>::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<ClusterType, Enable>::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<ClusterType, Enable>::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;
 }
 

include/aare/ClusterVector.hpp

@@ -32,7 +32,6 @@ class ClusterVector; // Forward declaration
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
           typename CoordType>
 class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
-    using value_type = T;
 
     std::byte *m_data{};
     size_t m_size{0};
@@ -49,6 +48,7 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
     constexpr static char m_fmt_base[] = "=h:x:\nh:y:\n({},{}){}:data:";
 
   public:
+    using value_type = T;
     using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
 
     /**
@@ -237,28 +237,6 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
         m_size = new_size;
     }
 
-    // TODO: Generalize !!!! Maybe move somewhere else
-    void apply_gain_map(const NDView<double> gain_map) {
-        // in principle we need to know the size of the image for this lookup
-        // TODO! check orientations
-        std::array<int64_t, 9> xcorr = {-1, 0, 1, -1, 0, 1, -1, 0, 1};
-        std::array<int64_t, 9> 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<T>(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;

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 <memory>
+
+namespace aare {
+
+class GainMap {
+
+  public:
+    explicit GainMap(const NDArray<double, 2> &gain_map)
+        : m_gain_map(gain_map) {};
+
+    explicit GainMap(const NDView<double, 2> gain_map) {
+        m_gain_map = NDArray<double, 2>(gain_map);
+    }
+
+    template <typename ClusterType,
+              typename = std::enable_if_t<is_cluster_v<ClusterType>>>
+    void apply_gain_map(ClusterVector<ClusterType> &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<ClusterType>::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<T>(cl.data[j] * m_gain_map(y, x));
+                }
+            } else {
+                memset(cl.data, 0,
+                       ClusterSizeX * ClusterSizeY *
+                           sizeof(T)); // clear edge clusters
+            }
+        }
+    }
+
+  private:
+    NDArray<double, 2> m_gain_map{};
+};
+
+} // end of namespace aare

src/ClusterVector.test.cpp

@@ -1,6 +1,7 @@
 #include "aare/ClusterVector.hpp"
 #include <cstdint>
 
+#include <catch2/catch_all.hpp>
 #include <catch2/catch_test_macros.hpp>
 #include <catch2/matchers/catch_matchers_floating_point.hpp>
 
@@ -184,3 +185,48 @@ TEST_CASE("Concatenate two cluster vectors where we need to allocate",
     REQUIRE(ptr[3].x == 16);
     REQUIRE(ptr[3].y == 17);
 }
+
+struct ClusterTestData {
+    int8_t ClusterSizeX;
+    int8_t ClusterSizeY;
+    std::vector<int64_t> index_map_x;
+    std::vector<int64_t> 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<int64_t> index_map_x(ClusterSizeX * ClusterSizeY);
+    std::vector<int64_t> 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);
+}