cpp Cluster and ClusterVector and ClusterFile are templated now, they support generic cluster types

include/aare/Cluster.hpp

@@ -0,0 +1,38 @@
+
+/************************************************
+ * @file Cluster.hpp
+ * @short definition of cluster, where CoordType (x,y) give
+ * the cluster center coordinates and data the actual cluster data
+ * cluster size is given as template parameters
+ ***********************************************/
+
+#pragma once
+
+#include <cstdint>
+#include <type_traits>
+
+namespace aare {
+
+// requires clause c++20 maybe update
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = int16_t,
+          typename Enable = std::enable_if_t<std::is_arithmetic_v<T> &&
+                                             std::is_integral_v<CoordType>>>
+struct Cluster {
+    CoordType x;
+    CoordType y;
+    T data[ClusterSizeX * ClusterSizeY];
+};
+
+// Type Traits for is_cluster_type
+template <typename T>
+struct is_cluster : std::false_type {}; // Default case: Not a Cluster
+
+// TODO: Do i need the require clause here as well?
+template <typename T, uint8_t X, uint8_t Y, typename CoordType>
+struct is_cluster<Cluster<T, X, Y, CoordType>> : std::true_type {}; // Cluster
+
+// helper
+template <typename T> constexpr bool is_cluster_v = is_cluster<T>::value;
+
+} // namespace aare

include/aare/ClusterFile.hpp

@@ -1,5 +1,6 @@
 #pragma once
 
+#include "aare/Cluster.hpp"
 #include "aare/ClusterVector.hpp"
 #include "aare/NDArray.hpp"
 #include "aare/defs.hpp"
@@ -8,14 +9,6 @@
 
 namespace aare {
 
-template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType = int16_t>
-struct Cluster {
-    CoordType x;
-    CoordType y;
-    T data[ClusterSizeX * ClusterSizeY];
-};
-
 typedef enum {
     cBottomLeft = 0,
     cBottomRight = 1,
@@ -59,6 +52,8 @@ uint32_t number_of_clusters
 ....
 */
 
+// TODO: change to support any type of clusters, e.g. header line with
+// clsuter_size_x, cluster_size_y,
 /**
  * @brief Class to read and write cluster files
  * Expects data to be laid out as:
@@ -71,6 +66,8 @@ uint32_t number_of_clusters
  *       uint32_t number_of_clusters
  *       etc.
  */
+template <typename ClusterType,
+          typename Enable = std::enable_if_t<is_cluster_v<ClusterType>, bool>>
 class ClusterFile {
     FILE *fp{};
     uint32_t m_num_left{};
@@ -97,9 +94,9 @@ class ClusterFile {
      * If EOF is reached the returned vector will have less than n_clusters
      * clusters
      */
-    ClusterVector<int32_t> read_clusters(size_t n_clusters);
+    ClusterVector<ClusterType> read_clusters(size_t n_clusters);
 
-    ClusterVector<int32_t> read_clusters(size_t n_clusters, ROI roi);
+    ClusterVector<ClusterType> read_clusters(size_t n_clusters, ROI roi);
 
     /**
      * @brief Read a single frame from the file and return the clusters. The
@@ -107,9 +104,9 @@ class ClusterFile {
      * @throws std::runtime_error if the file is not opened for reading or the
      * file pointer not at the beginning of a frame
      */
-    ClusterVector<int32_t> read_frame();
+    ClusterVector<ClusterType> read_frame();
 
-    void write_frame(const ClusterVector<int32_t> &clusters);
+    void write_frame(const ClusterVector<ClusterType> &clusters);
 
     // Need to be migrated to support NDArray and return a ClusterVector
     // std::vector<Cluster3x3>
@@ -130,20 +127,484 @@ class ClusterFile {
 
 int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad,
                  double *eta2x, double *eta2y, double *eta3x, double *eta3y);
-int analyze_cluster(Cluster3x3 &cl, int32_t *t2, int32_t *t3, char *quad,
+int analyze_cluster(Cluster<int32_t, 3, 3> &cl, int32_t *t2, int32_t *t3,
-                    double *eta2x, double *eta2y, double *eta3x, double *eta3y);
+                    char *quad, double *eta2x, double *eta2y, double *eta3x,
+                    double *eta3y);
 
-template <typename ClusterType>
+template <typename ClusterType,
+          typename = std::enable_if_t<is_cluster_v<ClusterType>>>
 NDArray<double, 2> calculate_eta2(ClusterVector<ClusterType> &clusters);
 
-template <typename T> Eta2 calculate_eta2(Cluster<T, 3, 3> &cl);
+// TODO: do we need rquire clauses?
+template <typename T> Eta2 calculate_eta2(const Cluster<T, 3, 3> &cl);
 
-Eta2 calculate_eta2(Cluster2x2 &cl);
+template <typename T> Eta2 calculate_eta2(const Cluster<T, 2, 2> &cl);
 
-template <typename ClusterType> Eta2 calculate_eta2(ClusterType &cl);
+template <typename ClusterType, std::enable_if_t<is_cluster_v<ClusterType>>>
+Eta2 calculate_eta2(const ClusterType &cl);
 
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
           typename CoordType>
-Eta2 calculate_eta2(Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl);
+Eta2 calculate_eta2(
+    const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl);
+
+template <typename ClusterType, typename Enable>
+ClusterFile<ClusterType, Enable>::ClusterFile(
+    const std::filesystem::path &fname, size_t chunk_size,
+    const std::string &mode)
+    : m_chunk_size(chunk_size), m_mode(mode) {
+
+    if (mode == "r") {
+        fp = fopen(fname.c_str(), "rb");
+        if (!fp) {
+            throw std::runtime_error("Could not open file for reading: " +
+                                     fname.string());
+        }
+    } else if (mode == "w") {
+        fp = fopen(fname.c_str(), "wb");
+        if (!fp) {
+            throw std::runtime_error("Could not open file for writing: " +
+                                     fname.string());
+        }
+    } else if (mode == "a") {
+        fp = fopen(fname.c_str(), "ab");
+        if (!fp) {
+            throw std::runtime_error("Could not open file for appending: " +
+                                     fname.string());
+        }
+    } else {
+        throw std::runtime_error("Unsupported mode: " + mode);
+    }
+}
+
+template <typename ClusterType, typename Enable>
+ClusterFile<ClusterType, Enable>::~ClusterFile() {
+    close();
+}
+
+template <typename ClusterType, typename Enable>
+void ClusterFile<ClusterType, Enable>::close() {
+    if (fp) {
+        fclose(fp);
+        fp = nullptr;
+    }
+}
+
+// TODO generally supported for all clsuter types
+template <typename ClusterType, typename Enable>
+void ClusterFile<ClusterType, Enable>::write_frame(
+    const ClusterVector<ClusterType> &clusters) {
+    if (m_mode != "w" && m_mode != "a") {
+        throw std::runtime_error("File not opened for writing");
+    }
+    if (!(clusters.cluster_size_x() == 3) &&
+        !(clusters.cluster_size_y() == 3)) {
+        throw std::runtime_error("Only 3x3 clusters are supported");
+    }
+    int32_t frame_number = clusters.frame_number();
+    fwrite(&frame_number, sizeof(frame_number), 1, fp);
+    uint32_t n_clusters = clusters.size();
+    fwrite(&n_clusters, sizeof(n_clusters), 1, fp);
+    fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp);
+}
+
+template <typename ClusterType, typename Enable>
+ClusterVector<ClusterType>
+ClusterFile<ClusterType, Enable>::read_clusters(size_t n_clusters) {
+    if (m_mode != "r") {
+        throw std::runtime_error("File not opened for reading");
+    }
+
+    ClusterVector<ClusterType> clusters(n_clusters);
+
+    int32_t iframe = 0; // frame number needs to be 4 bytes!
+    size_t nph_read = 0;
+    uint32_t nn = m_num_left;
+    uint32_t nph = m_num_left; // number of clusters in frame needs to be 4
+
+    // auto buf = reinterpret_cast<Cluster3x3 *>(clusters.data());
+    auto buf = clusters.data();
+    // if there are photons left from previous frame read them first
+    if (nph) {
+        if (nph > n_clusters) {
+            // if we have more photons left in the frame then photons to read we
+            // read directly the requested number
+            nn = n_clusters;
+        } else {
+            nn = nph;
+        }
+        nph_read += fread((buf + nph_read * clusters.item_size()),
+                          clusters.item_size(), nn, fp);
+        m_num_left = nph - nn; // write back the number of photons left
+    }
+
+    if (nph_read < n_clusters) {
+        // keep on reading frames and photons until reaching n_clusters
+        while (fread(&iframe, sizeof(iframe), 1, fp)) {
+            // read number of clusters in frame
+            if (fread(&nph, sizeof(nph), 1, fp)) {
+                if (nph > (n_clusters - nph_read))
+                    nn = n_clusters - nph_read;
+                else
+                    nn = nph;
+
+                nph_read += fread((buf + nph_read * clusters.item_size()),
+                                  clusters.item_size(), nn, fp);
+                m_num_left = nph - nn;
+            }
+            if (nph_read >= n_clusters)
+                break;
+        }
+    }
+
+    // Resize the vector to the number of clusters.
+    // No new allocation, only change bounds.
+    clusters.resize(nph_read);
+    return clusters;
+}
+
+template <typename ClusterType, typename Enable>
+ClusterVector<ClusterType>
+ClusterFile<ClusterType, Enable>::read_clusters(size_t n_clusters, ROI roi) {
+    if (m_mode != "r") {
+        throw std::runtime_error("File not opened for reading");
+    }
+
+    ClusterVector<ClusterType> clusters;
+    clusters.reserve(n_clusters);
+
+    int32_t iframe = 0; // frame number needs to be 4 bytes!
+    size_t nph_read = 0;
+    uint32_t nn = m_num_left;
+    uint32_t nph = m_num_left; // number of clusters in frame needs to be 4
+
+    // auto buf = reinterpret_cast<Cluster3x3 *>(clusters.data());
+    // auto buf = clusters.data();
+
+    ClusterType tmp; // this would break if the cluster size changes
+
+    // if there are photons left from previous frame read them first
+    if (nph) {
+        if (nph > n_clusters) {
+            // if we have more photons left in the frame then photons to read we
+            // read directly the requested number
+            nn = n_clusters;
+        } else {
+            nn = nph;
+        }
+        // Read one cluster, in the ROI push back
+        //  nph_read += fread((buf + nph_read*clusters.item_size()),
+        //                    clusters.item_size(), nn, fp);
+        for (size_t i = 0; i < nn; i++) {
+            fread(&tmp, sizeof(tmp), 1, fp);
+            if (tmp.x >= roi.xmin && tmp.x <= roi.xmax && tmp.y >= roi.ymin &&
+                tmp.y <= roi.ymax) {
+                clusters.push_back(tmp.x, tmp.y,
+                                   reinterpret_cast<std::byte *>(tmp.data));
+                nph_read++;
+            }
+        }
+
+        m_num_left = nph - nn; // write back the number of photons left
+    }
+
+    if (nph_read < n_clusters) {
+        // keep on reading frames and photons until reaching n_clusters
+        while (fread(&iframe, sizeof(iframe), 1, fp)) {
+            // read number of clusters in frame
+            if (fread(&nph, sizeof(nph), 1, fp)) {
+                if (nph > (n_clusters - nph_read))
+                    nn = n_clusters - nph_read;
+                else
+                    nn = nph;
+
+                // nph_read += fread((buf + nph_read*clusters.item_size()),
+                //                   clusters.item_size(), nn, fp);
+                for (size_t i = 0; i < nn; i++) {
+                    fread(&tmp, sizeof(tmp), 1, fp);
+                    if (tmp.x >= roi.xmin && tmp.x <= roi.xmax &&
+                        tmp.y >= roi.ymin && tmp.y <= roi.ymax) {
+                        clusters.push_back(
+                            tmp.x, tmp.y,
+                            reinterpret_cast<std::byte *>(tmp.data));
+                        nph_read++;
+                    }
+                }
+                m_num_left = nph - nn;
+            }
+            if (nph_read >= n_clusters)
+                break;
+        }
+    }
+
+    // Resize the vector to the number of clusters.
+    // No new allocation, only change bounds.
+    clusters.resize(nph_read);
+    return clusters;
+}
+
+template <typename ClusterType, typename Enable>
+ClusterVector<ClusterType> ClusterFile<ClusterType, Enable>::read_frame() {
+    if (m_mode != "r") {
+        throw std::runtime_error("File not opened for reading");
+    }
+    if (m_num_left) {
+        throw std::runtime_error(
+            "There are still photons left in the last frame");
+    }
+    int32_t frame_number;
+    if (fread(&frame_number, sizeof(frame_number), 1, fp) != 1) {
+        throw std::runtime_error("Could not read frame number");
+    }
+
+    int32_t n_clusters; // Saved as 32bit integer in the cluster file
+    if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) {
+        throw std::runtime_error("Could not read number of clusters");
+    }
+    // std::vector<Cluster3x3> clusters(n_clusters);
+    ClusterVector<ClusterType> clusters(n_clusters);
+    clusters.set_frame_number(frame_number);
+
+    if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) !=
+        static_cast<size_t>(n_clusters)) {
+        throw std::runtime_error("Could not read clusters");
+    }
+    clusters.resize(n_clusters);
+    return clusters;
+}
+
+template <typename ClusterType, std::enable_if_t<is_cluster_v<ClusterType>>>
+NDArray<double, 2> calculate_eta2(const ClusterVector<ClusterType> &clusters) {
+    // TOTO! make work with 2x2 clusters
+    NDArray<double, 2> eta2({static_cast<int64_t>(clusters.size()), 2});
+
+    for (size_t i = 0; i < clusters.size(); i++) {
+        auto e = calculate_eta2<ClusterType>(clusters.at(i));
+        eta2(i, 0) = e.x;
+        eta2(i, 1) = e.y;
+    }
+
+    return eta2;
+}
+
+/**
+ * @brief Calculate the eta2 values for a generic sized cluster and return them
+ * in a Eta2 struct containing etay, etax and the index of the respective 2x2
+ * subcluster.
+ */
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType>
+Eta2 calculate_eta2(
+    const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
+    Eta2 eta{};
+
+    // TODO loads of overhead for a 2x2 clsuter maybe keep 2x2 calculation
+    constexpr size_t num_2x2_subclusters =
+        (ClusterSizeX - 1) * (ClusterSizeY - 1);
+    std::array<T, num_2x2_subclusters> sum_2x2_subcluster;
+    for (size_t i = 0; i < ClusterSizeY - 1; ++i) {
+        for (size_t j = 0; j < ClusterSizeX - 1; ++j)
+            sum_2x2_subcluster[i * (ClusterSizeX - 1) + j] =
+                cl.data[i * ClusterSizeX + j] +
+                cl.data[i * ClusterSizeX + j + 1] +
+                cl.data[(i + 1) * ClusterSizeX + j] +
+                cl.data[(i + 1) * ClusterSizeX + j + 1];
+    }
+
+    auto c =
+        std::max_element(sum_2x2_subcluster.begin(), sum_2x2_subcluster.end()) -
+        sum_2x2_subcluster.begin();
+
+    eta.sum = sum_2x2_subcluster[c];
+
+    eta.x = static_cast<double>(cl.data[(c + 1) * ClusterSizeX + 1]) /
+            (cl.data[0] + cl.data[1]);
+
+    size_t index_top_left_2x2_subcluster =
+        (int(c / (ClusterSizeX - 1)) + 1) * ClusterSizeX +
+        c % (ClusterSizeX - 1) * 2 + 1;
+    if ((cl.data[index_top_left_2x2_subcluster] +
+         cl.data[index_top_left_2x2_subcluster - 1]) != 0)
+        eta.x =
+            static_cast<double>(cl.data[index_top_left_2x2_subcluster] /
+                                (cl.data[index_top_left_2x2_subcluster] +
+                                 cl.data[index_top_left_2x2_subcluster - 1]));
+
+    if ((cl.data[index_top_left_2x2_subcluster] +
+         cl.data[index_top_left_2x2_subcluster - ClusterSizeX]) != 0)
+        eta.y = static_cast<double>(
+            cl.data[index_top_left_2x2_subcluster] /
+            (cl.data[index_top_left_2x2_subcluster] +
+             cl.data[index_top_left_2x2_subcluster - ClusterSizeX]));
+
+    eta.c = c; // TODO only supported for 2x2 and 3x3 clusters -> at least no
+               // underyling enum class
+    return eta;
+}
+
+/**
+ * @brief Calculate the eta2 values for a 3x3 cluster and return them in a Eta2
+ * struct containing etay, etax and the corner of the cluster.
+ */
+template <typename T> Eta2 calculate_eta2(const Cluster<T, 3, 3> &cl) {
+    Eta2 eta{};
+
+    std::array<T, 4> tot2;
+    tot2[0] = cl.data[0] + cl.data[1] + cl.data[3] + cl.data[4];
+    tot2[1] = cl.data[1] + cl.data[2] + cl.data[4] + cl.data[5];
+    tot2[2] = cl.data[3] + cl.data[4] + cl.data[6] + cl.data[7];
+    tot2[3] = cl.data[4] + cl.data[5] + cl.data[7] + cl.data[8];
+
+    auto c = std::max_element(tot2.begin(), tot2.end()) - tot2.begin();
+    eta.sum = tot2[c];
+    switch (c) {
+    case cBottomLeft:
+        if ((cl.data[3] + cl.data[4]) != 0)
+            eta.x = static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
+        if ((cl.data[1] + cl.data[4]) != 0)
+            eta.y = static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
+        eta.c = cBottomLeft;
+        break;
+    case cBottomRight:
+        if ((cl.data[2] + cl.data[5]) != 0)
+            eta.x = static_cast<double>(cl.data[5]) / (cl.data[4] + cl.data[5]);
+        if ((cl.data[1] + cl.data[4]) != 0)
+            eta.y = static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
+        eta.c = cBottomRight;
+        break;
+    case cTopLeft:
+        if ((cl.data[7] + cl.data[4]) != 0)
+            eta.x = static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
+        if ((cl.data[7] + cl.data[4]) != 0)
+            eta.y = static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
+        eta.c = cTopLeft;
+        break;
+    case cTopRight:
+        if ((cl.data[5] + cl.data[4]) != 0)
+            eta.x = static_cast<double>(cl.data[5]) / (cl.data[5] + cl.data[4]);
+        if ((cl.data[7] + cl.data[4]) != 0)
+            eta.y = static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
+        eta.c = cTopRight;
+        break;
+    }
+    return eta;
+}
+
+template <typename T> Eta2 calculate_eta2(const Cluster<T, 2, 2> &cl) {
+    Eta2 eta{};
+
+    eta.x = static_cast<double>(cl.data[1]) / (cl.data[0] + cl.data[1]);
+    eta.y = static_cast<double>(cl.data[2]) / (cl.data[0] + cl.data[2]);
+    eta.sum = cl.data[0] + cl.data[1] + cl.data[2] + cl.data[3];
+    eta.c = cBottomLeft; // TODO! This is not correct, but need to put something
+    return eta;
+}
+
+// TODO complicated API simplify?
+int analyze_cluster(Cluster<int32_t, 3, 3> &cl, int32_t *t2, int32_t *t3,
+                    char *quad, double *eta2x, double *eta2y, double *eta3x,
+                    double *eta3y) {
+
+    return analyze_data(cl.data, t2, t3, quad, eta2x, eta2y, eta3x, eta3y);
+}
+
+int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad,
+                 double *eta2x, double *eta2y, double *eta3x, double *eta3y) {
+
+    int ok = 1;
+
+    int32_t tot2[4];
+    int32_t t2max = 0;
+    char c = 0;
+    int32_t val, tot3;
+
+    tot3 = 0;
+    for (int i = 0; i < 4; i++)
+        tot2[i] = 0;
+
+    for (int ix = 0; ix < 3; ix++) {
+        for (int iy = 0; iy < 3; iy++) {
+            val = data[iy * 3 + ix];
+            //	printf ("%d ",data[iy * 3 + ix]);
+            tot3 += val;
+            if (ix <= 1 && iy <= 1)
+                tot2[cBottomLeft] += val;
+            if (ix >= 1 && iy <= 1)
+                tot2[cBottomRight] += val;
+            if (ix <= 1 && iy >= 1)
+                tot2[cTopLeft] += val;
+            if (ix >= 1 && iy >= 1)
+                tot2[cTopRight] += val;
+        }
+        //	printf ("\n");
+    }
+    // printf ("\n");
+
+    if (t2 || quad) {
+
+        t2max = tot2[0];
+        c = cBottomLeft;
+        for (int i = 1; i < 4; i++) {
+            if (tot2[i] > t2max) {
+                t2max = tot2[i];
+                c = i;
+            }
+        }
+        // printf("*** %d %d %d %d --
+        // %d\n",tot2[0],tot2[1],tot2[2],tot2[3],t2max);
+        if (quad)
+            *quad = c;
+        if (t2)
+            *t2 = t2max;
+    }
+
+    if (t3)
+        *t3 = tot3;
+
+    if (eta2x || eta2y) {
+        if (eta2x)
+            *eta2x = 0;
+        if (eta2y)
+            *eta2y = 0;
+        switch (c) {
+        case cBottomLeft:
+            if (eta2x && (data[3] + data[4]) != 0)
+                *eta2x = static_cast<double>(data[4]) / (data[3] + data[4]);
+            if (eta2y && (data[1] + data[4]) != 0)
+                *eta2y = static_cast<double>(data[4]) / (data[1] + data[4]);
+            break;
+        case cBottomRight:
+            if (eta2x && (data[2] + data[5]) != 0)
+                *eta2x = static_cast<double>(data[5]) / (data[4] + data[5]);
+            if (eta2y && (data[1] + data[4]) != 0)
+                *eta2y = static_cast<double>(data[4]) / (data[1] + data[4]);
+            break;
+        case cTopLeft:
+            if (eta2x && (data[7] + data[4]) != 0)
+                *eta2x = static_cast<double>(data[4]) / (data[3] + data[4]);
+            if (eta2y && (data[7] + data[4]) != 0)
+                *eta2y = static_cast<double>(data[7]) / (data[7] + data[4]);
+            break;
+        case cTopRight:
+            if (eta2x && t2max != 0)
+                *eta2x = static_cast<double>(data[5]) / (data[5] + data[4]);
+            if (eta2y && t2max != 0)
+                *eta2y = static_cast<double>(data[7]) / (data[7] + data[4]);
+            break;
+        default:;
+        }
+    }
+
+    if (eta3x || eta3y) {
+        if (eta3x && (data[3] + data[4] + data[5]) != 0)
+            *eta3x = static_cast<double>(-data[3] + data[3 + 2]) /
+                     (data[3] + data[4] + data[5]);
+        if (eta3y && (data[1] + data[4] + data[7]) != 0)
+            *eta3y = static_cast<double>(-data[1] + data[2 * 3 + 1]) /
+                     (data[1] + data[4] + data[7]);
+    }
+
+    return ok;
+}
 
 } // namespace aare

include/aare/ClusterVector.hpp

@@ -1,4 +1,5 @@
 #pragma once
+#include "aare/Cluster.hpp" //TODO maybe store in seperate file !!!
 #include <algorithm>
 #include <array>
 #include <cstddef>
@@ -10,7 +11,9 @@
 
 namespace aare {
 
-template <typename ClusterType> class ClusterVector; // Forward declaration
+template <typename ClusterType,
+          typename = std::enable_if_t<is_cluster_v<ClusterType>>>
+class ClusterVector; // Forward declaration
 
 /**
  * @brief ClusterVector is a container for clusters of various sizes. It uses a
@@ -44,12 +47,10 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
     constexpr static char m_fmt_base[] = "=h:x:\nh:y:\n({},{}){}:data:";
 
   public:
-    using ClusterType = Cluster<T, SizeX, SizeY>;
+    using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
 
     /**
      * @brief Construct a new ClusterVector object
-     * @param cluster_size_x size of the cluster in x direction
-     * @param cluster_size_y size of the cluster in y direction
      * @param capacity initial capacity of the buffer in number of clusters
      * @param frame_number frame number of the clusters. Default is 0, which is
      * also used to indicate that the clusters come from many frames
@@ -184,6 +185,10 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
      */
     size_t size() const { return m_size; }
 
+    uint8_t cluster_size_x() const { return ClusterSizeX; }
+
+    uint8_t cluster_size_y() const { return ClusterSizeY; }
+
     /**
      * @brief Return the capacity of the buffer in number of clusters. This is
      * the number of clusters that can be stored in the current buffer without

include/aare/Interpolator.hpp

@@ -1,29 +1,35 @@
 #pragma once
+
+#include "aare/Cluster.hpp"
+#include "aare/ClusterFile.hpp" //Cluster_3x3
+#include "aare/ClusterVector.hpp"
 #include "aare/NDArray.hpp"
 #include "aare/NDView.hpp"
-#include "aare/ClusterVector.hpp"
+namespace aare {
-#include "aare/ClusterFile.hpp" //Cluster_3x3
-namespace aare{
 
-struct Photon{
+struct Photon {
     double x;
     double y;
     double energy;
 };
 
-class Interpolator{
+class Interpolator {
     NDArray<double, 3> m_ietax;
     NDArray<double, 3> m_ietay;
 
     NDArray<double, 1> m_etabinsx;
     NDArray<double, 1> m_etabinsy;
     NDArray<double, 1> m_energy_bins;
-    public:
-        Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins, NDView<double, 1> ybins, NDView<double, 1> ebins);
-        NDArray<double, 3> get_ietax(){return m_ietax;}
-        NDArray<double, 3> get_ietay(){return m_ietay;}
 
-        std::vector<Photon> interpolate(const ClusterVector<int32_t>& clusters);
+  public:
+    Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
+                 NDView<double, 1> ybins, NDView<double, 1> ebins);
+    NDArray<double, 3> get_ietax() { return m_ietax; }
+    NDArray<double, 3> get_ietay() { return m_ietay; }
+
+    template <typename ClusterType,
+              typename = std::enable_if_t<is_cluster_v<ClusterType>>>
+    std::vector<Photon> interpolate(const ClusterVector<ClusterType> &clusters);
 };
 
 } // namespace aare

src/ClusterFile.cpp

@@ -4,8 +4,11 @@
 
 namespace aare {
 
-ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size,
+template <typename ClusterType,
-                         const std::string &mode)
+          typename = std::enable_if_t<is_cluster_v<ClusterType>>>
+ClusterFile<ClusterType>::ClusterFile(const std::filesystem::path &fname,
+                                      size_t chunk_size,
+                                      const std::string &mode)
     : m_chunk_size(chunk_size), m_mode(mode) {
 
     if (mode == "r") {
@@ -31,16 +34,21 @@ ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size,
     }
 }
 
-ClusterFile::~ClusterFile() { close(); }
+template <typename ClusterType> ClusterFile<ClusterType>::~ClusterFile() {
+    close();
+}
 
-void ClusterFile::close() {
+template <typename ClusterType> void ClusterFile<ClusterType>::close() {
     if (fp) {
         fclose(fp);
         fp = nullptr;
     }
 }
 
-void ClusterFile::write_frame(const ClusterVector<int32_t> &clusters) {
+// TODO generally supported for all clsuter types
+template <typename ClusterType>
+void ClusterFile<ClusterType>::write_frame(
+    const ClusterVector<ClusterType> &clusters) {
     if (m_mode != "w" && m_mode != "a") {
         throw std::runtime_error("File not opened for writing");
     }
@@ -55,12 +63,14 @@ void ClusterFile::write_frame(const ClusterVector<int32_t> &clusters) {
     fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp);
 }
 
-ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters) {
+template <typename ClusterType>
+ClusterVector<ClusterType>
+ClusterFile<ClusterType>::read_clusters(size_t n_clusters) {
     if (m_mode != "r") {
         throw std::runtime_error("File not opened for reading");
     }
 
-    ClusterVector<int32_t> clusters(3, 3, n_clusters);
+    ClusterVector<ClusterType> clusters(n_clusters);
 
     int32_t iframe = 0; // frame number needs to be 4 bytes!
     size_t nph_read = 0;
@@ -108,12 +118,14 @@ ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters) {
     return clusters;
 }
 
-ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters, ROI roi) {
+template <typename ClusterType>
+ClusterVector<ClusterType>
+ClusterFile<ClusterType>::read_clusters(size_t n_clusters, ROI roi) {
     if (m_mode != "r") {
         throw std::runtime_error("File not opened for reading");
     }
 
-    ClusterVector<int32_t> clusters(3, 3);
+    ClusterVector<ClusterType> clusters;
     clusters.reserve(n_clusters);
 
     int32_t iframe = 0; // frame number needs to be 4 bytes!
@@ -124,7 +136,7 @@ ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters, ROI roi) {
     // auto buf = reinterpret_cast<Cluster3x3 *>(clusters.data());
     // auto buf = clusters.data();
 
-    Cluster3x3 tmp; // this would break if the cluster size changes
+    ClusterType tmp; // this would break if the cluster size changes
 
     // if there are photons left from previous frame read them first
     if (nph) {
@@ -186,7 +198,8 @@ ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters, ROI roi) {
     return clusters;
 }
 
-ClusterVector<int32_t> ClusterFile::read_frame() {
+template <typename ClusterType>
+ClusterVector<ClusterType> ClusterFile<ClusterType>::read_frame() {
     if (m_mode != "r") {
         throw std::runtime_error("File not opened for reading");
     }
@@ -204,7 +217,7 @@ ClusterVector<int32_t> ClusterFile::read_frame() {
         throw std::runtime_error("Could not read number of clusters");
     }
     // std::vector<Cluster3x3> clusters(n_clusters);
-    ClusterVector<int32_t> clusters(3, 3, n_clusters);
+    ClusterVector<ClusterType> clusters(n_clusters);
     clusters.set_frame_number(frame_number);
 
     if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) !=
@@ -372,8 +385,9 @@ Eta2 calculate_eta2(Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
     Eta2 eta{};
 
     // TODO loads of overhead for a 2x2 clsuter maybe keep 2x2 calculation
-    size_t num_2x2_subclusters = (ClusterSizeX - 1) * (ClusterSizeY - 1);
+    constexpr size_t num_2x2_subclusters =
-    std::array<int32_t, num_2x2_subclusters> sum_2x2_subcluster;
+        (ClusterSizeX - 1) * (ClusterSizeY - 1);
+    std::array<T, num_2x2_subclusters> sum_2x2_subcluster;
     for (size_t i = 0; i < ClusterSizeY - 1; ++i) {
         for (size_t j = 0; j < ClusterSizeX - 1; ++j)
             sum_2x2_subcluster[i * (ClusterSizeX - 1) + j] =
@@ -383,9 +397,9 @@ Eta2 calculate_eta2(Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
                 cl.data[(i + 1) * ClusterSizeX + j + 1];
     }
 
-    auto c = std::max_element(sum_2x2_subclusters.begin(),
+    auto c =
-                              sum_2x2_subcluster.end()) -
+        std::max_element(sum_2x2_subcluster.begin(), sum_2x2_subcluster.end()) -
-             sum_2x2_subcluster.begin();
+        sum_2x2_subcluster.begin();
 
     eta.sum = sum_2x2_subcluster[c];
 
@@ -458,7 +472,6 @@ template <typename T> Eta2 calculate_eta2(Cluster<T, 3, 3> &cl) {
             eta.y = static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
         eta.c = cTopRight;
         break;
-        // no default to allow compiler to warn about missing cases
     }
     return eta;
 }
@@ -473,8 +486,9 @@ template <typename T> Eta2 calculate_eta2(Cluster<T, 2, 2> &cl) {
     return eta;
 }
 
-int analyze_cluster(Cluster3x3 &cl, int32_t *t2, int32_t *t3, char *quad,
+// TODO complicated API simplify?
-                    double *eta2x, double *eta2y, double *eta3x,
+int analyze_cluster(Cluster<int32_t, 3, 3> &cl, int32_t *t2, int32_t *t3,
+                    char *quad, double *eta2x, double *eta2y, double *eta3x,
                     double *eta3y) {
 
     return analyze_data(cl.data, t2, t3, quad, eta2x, eta2y, eta3x, eta3y);

src/Interpolator.cpp

@@ -5,7 +5,8 @@ namespace aare {
 
 Interpolator::Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
                            NDView<double, 1> ybins, NDView<double, 1> ebins)
-    : m_ietax(etacube), m_ietay(etacube), m_etabinsx(xbins), m_etabinsy(ybins), m_energy_bins(ebins) {
+    : m_ietax(etacube), m_ietay(etacube), m_etabinsx(xbins), m_etabinsy(ybins),
+      m_energy_bins(ebins) {
     if (etacube.shape(0) != xbins.size() || etacube.shape(1) != ybins.size() ||
         etacube.shape(2) != ebins.size()) {
         throw std::invalid_argument(
@@ -51,35 +52,37 @@ Interpolator::Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
             }
         }
     }
-
 }
 
-std::vector<Photon> Interpolator::interpolate(const ClusterVector<int32_t>& clusters) {
+template <typename ClusterType,
+          typename = std::enable_if_t<is_cluster_v<ClusterType>>>
+std::vector<Photon>
+Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) {
     std::vector<Photon> photons;
     photons.reserve(clusters.size());
 
     if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) {
-        for (size_t i = 0; i<clusters.size(); i++){
+        for (size_t i = 0; i < clusters.size(); i++) {
+
+            auto cluster = clusters.at(i);
+            Eta2 eta = calculate_eta2(cluster);
 
-            auto cluster = clusters.at<Cluster3x3>(i);
-            Eta2 eta= calculate_eta2(cluster);
-    
             Photon photon;
             photon.x = cluster.x;
             photon.y = cluster.y;
             photon.energy = eta.sum;
-    
+
             // auto ie = nearest_index(m_energy_bins, photon.energy)-1;
             // auto ix = nearest_index(m_etabinsx, eta.x)-1;
-            // auto iy = nearest_index(m_etabinsy, eta.y)-1;   
+            // auto iy = nearest_index(m_etabinsy, eta.y)-1;
-            //Finding the index of the last element that is smaller
+            // Finding the index of the last element that is smaller
-            //should work fine as long as we have many bins
+            // should work fine as long as we have many bins
             auto ie = last_smaller(m_energy_bins, photon.energy);
             auto ix = last_smaller(m_etabinsx, eta.x);
-            auto iy = last_smaller(m_etabinsy, eta.y); 
+            auto iy = last_smaller(m_etabinsy, eta.y);
 
             // fmt::print("ex: {}, ix: {}, iy: {}\n", ie, ix, iy);
-            
+
             double dX, dY;
             int ex, ey;
             // cBottomLeft = 0,
@@ -100,44 +103,47 @@ std::vector<Photon> Interpolator::interpolate(const ClusterVector<int32_t>& clus
                 dY = -1.;
                 break;
             case cBottomRight:
-                dX = 0.; 
+                dX = 0.;
                 dY = -1.;
                 break;
             }
-            photon.x += m_ietax(ix, iy, 0)*2 + dX;
+            photon.x += m_ietax(ix, iy, 0) * 2 + dX;
-            photon.y += m_ietay(ix, iy, 0)*2 + dY;
+            photon.y += m_ietay(ix, iy, 0) * 2 + dY;
             photons.push_back(photon);
         }
-    }else if(clusters.cluster_size_x() == 2 || clusters.cluster_size_y() == 2){
+    } else if (clusters.cluster_size_x() == 2 ||
-        for (size_t i = 0; i<clusters.size(); i++){
+               clusters.cluster_size_y() == 2) {
-            auto cluster = clusters.at<Cluster2x2>(i);
+        for (size_t i = 0; i < clusters.size(); i++) {
-            Eta2 eta= calculate_eta2(cluster);
+            auto cluster = clusters.at(i);
-    
+            Eta2 eta = calculate_eta2(cluster);
+
             Photon photon;
             photon.x = cluster.x;
             photon.y = cluster.y;
             photon.energy = eta.sum;
-    
+
-            //Now do some actual interpolation. 
+            // Now do some actual interpolation.
-            //Find which energy bin the cluster is in
+            // Find which energy bin the cluster is in
-            // auto ie = nearest_index(m_energy_bins, photon.energy)-1;
+            //  auto ie = nearest_index(m_energy_bins, photon.energy)-1;
-            // auto ix = nearest_index(m_etabinsx, eta.x)-1;
+            //  auto ix = nearest_index(m_etabinsx, eta.x)-1;
-            // auto iy = nearest_index(m_etabinsy, eta.y)-1;  
+            //  auto iy = nearest_index(m_etabinsy, eta.y)-1;
-            //Finding the index of the last element that is smaller
+            // Finding the index of the last element that is smaller
-            //should work fine as long as we have many bins
+            // should work fine as long as we have many bins
             auto ie = last_smaller(m_energy_bins, photon.energy);
             auto ix = last_smaller(m_etabinsx, eta.x);
-            auto iy = last_smaller(m_etabinsy, eta.y); 
+            auto iy = last_smaller(m_etabinsy, eta.y);
 
-            photon.x += m_ietax(ix, iy, 0)*2; //eta goes between 0 and 1 but we could move the hit anywhere in the 2x2
+            photon.x +=
-            photon.y += m_ietay(ix, iy, 0)*2;
+                m_ietax(ix, iy, 0) * 2; // eta goes between 0 and 1 but we could
+                                        // move the hit anywhere in the 2x2
+            photon.y += m_ietay(ix, iy, 0) * 2;
             photons.push_back(photon);
         }
-        
+
-    }else{
+    } else {
-        throw std::runtime_error("Only 3x3 and 2x2 clusters are supported for interpolation");
+        throw std::runtime_error(
+            "Only 3x3 and 2x2 clusters are supported for interpolation");
     }
-  
 
     return photons;
 }