removed extra sum after merge

- Change of pointer size caused out of bounds write
- UB to write to memory reserved by std::vector::reserver --> allocate
dummy clusters by using resize instead
   - but now we can't reserve like we want to, need a fix. 
- format string not working, fixed

CMakeLists.txt

@@ -384,7 +384,6 @@ set(SourceFiles
     ${CMAKE_CURRENT_SOURCE_DIR}/src/utils/task.cpp
 )
 
-
 add_library(aare_core STATIC ${SourceFiles})
 target_include_directories(aare_core PUBLIC 
     "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>" 

include/aare/CalculateEta.hpp

@@ -33,7 +33,7 @@ template <typename T> struct Eta2 {
 };
 
 /**
- * @brief Calculate the eta2 values for all clusters in a Clsutervector
+ * @brief Calculate the eta2 values for all clusters in a Clustervector
  */
 template <typename ClusterType,
           typename = std::enable_if_t<is_cluster_v<ClusterType>>>

include/aare/ClusterFile.hpp

@@ -39,9 +39,10 @@ template <typename ClusterType,
           typename Enable = std::enable_if_t<is_cluster_v<ClusterType>>>
 class ClusterFile {
     FILE *fp{};
+    const std::string m_filename{};
     uint32_t m_num_left{};    /*Number of photons left in frame*/
     size_t m_chunk_size{};    /*Number of clusters to read at a time*/
-    const std::string m_mode; /*Mode to open the file in*/
+    std::string m_mode;       /*Mode to open the file in*/
     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*/
@@ -115,6 +116,11 @@ class ClusterFile {
      */
     void close();
 
+    /** @brief Open the file in specific mode
+     *
+     */
+    void open(const std::string &mode);
+
   private:
     ClusterVector<ClusterType> read_clusters_with_cut(size_t n_clusters);
     ClusterVector<ClusterType> read_clusters_without_cut(size_t n_clusters);
@@ -128,25 +134,25 @@ 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) {
+    : m_filename(fname.string()), m_chunk_size(chunk_size), m_mode(mode) {
 
     if (mode == "r") {
-        fp = fopen(fname.c_str(), "rb");
+        fp = fopen(m_filename.c_str(), "rb");
         if (!fp) {
             throw std::runtime_error("Could not open file for reading: " +
-                                     fname.string());
+                                     m_filename);
         }
     } else if (mode == "w") {
-        fp = fopen(fname.c_str(), "wb");
+        fp = fopen(m_filename.c_str(), "wb");
         if (!fp) {
             throw std::runtime_error("Could not open file for writing: " +
-                                     fname.string());
+                                     m_filename);
         }
     } else if (mode == "a") {
-        fp = fopen(fname.c_str(), "ab");
+        fp = fopen(m_filename.c_str(), "ab");
         if (!fp) {
             throw std::runtime_error("Could not open file for appending: " +
-                                     fname.string());
+                                     m_filename);
         }
     } else {
         throw std::runtime_error("Unsupported mode: " + mode);
@@ -165,6 +171,39 @@ void ClusterFile<ClusterType, Enable>::close() {
         fp = nullptr;
     }
 }
+
+template <typename ClusterType, typename Enable>
+void ClusterFile<ClusterType, Enable>::open(const std::string &mode) {
+    if (fp) {
+        close();
+    }
+
+    if (mode == "r") {
+        fp = fopen(m_filename.c_str(), "rb");
+        if (!fp) {
+            throw std::runtime_error("Could not open file for reading: " +
+                                     m_filename);
+        }
+        m_mode = "r";
+    } else if (mode == "w") {
+        fp = fopen(m_filename.c_str(), "wb");
+        if (!fp) {
+            throw std::runtime_error("Could not open file for writing: " +
+                                     m_filename);
+        }
+        m_mode = "w";
+    } else if (mode == "a") {
+        fp = fopen(m_filename.c_str(), "ab");
+        if (!fp) {
+            throw std::runtime_error("Could not open file for appending: " +
+                                     m_filename);
+        }
+        m_mode = "a";
+    } else {
+        throw std::runtime_error("Unsupported mode: " + mode);
+    }
+}
+
 template <typename ClusterType, typename Enable>
 void ClusterFile<ClusterType, Enable>::set_roi(ROI roi) {
     m_roi = roi;
@@ -197,10 +236,7 @@ void ClusterFile<ClusterType, Enable>::write_frame(
     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();
@@ -270,7 +306,7 @@ ClusterFile<ClusterType, Enable>::read_clusters_without_cut(size_t n_clusters) {
         }
     }
 
-    // Resize the vector to the number of clusters.
+    // Resize the vector to the number o f clusters.
     // No new allocation, only change bounds.
     clusters.resize(nph_read);
     if (m_gain_map)
@@ -282,7 +318,7 @@ template <typename ClusterType, typename Enable>
 ClusterVector<ClusterType>
 ClusterFile<ClusterType, Enable>::read_clusters_with_cut(size_t n_clusters) {
     ClusterVector<ClusterType> clusters;
-    clusters.resize(n_clusters);
+    clusters.reserve(n_clusters);
 
     // if there are photons left from previous frame read them first
     if (m_num_left) {
@@ -375,11 +411,13 @@ ClusterFile<ClusterType, Enable>::read_frame_without_cut() {
     ClusterVector<ClusterType> clusters(n_clusters);
     clusters.set_frame_number(frame_number);
 
+    clusters.resize(n_clusters);
+
     if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) !=
         static_cast<size_t>(n_clusters)) {
         throw std::runtime_error(LOCATION + "Could not read clusters");
     }
-    clusters.resize(n_clusters);
+
     if (m_gain_map)
         m_gain_map->apply_gain_map(clusters);
     return clusters;

include/aare/ClusterFileV2.hpp

@@ -1,154 +0,0 @@
-#pragma once
-#include "aare/core/defs.hpp"
-#include <filesystem>
-#include <fmt/format.h>
-#include <string>
-
-namespace aare {
-struct ClusterHeader {
-    int32_t frame_number;
-    int32_t n_clusters;
-    std::string to_string() const {
-        return "frame_number: " + std::to_string(frame_number) +
-               ", n_clusters: " + std::to_string(n_clusters);
-    }
-};
-
-struct ClusterV2_ {
-    int16_t x;
-    int16_t y;
-    std::array<int32_t, 9> data;
-    std::string to_string(bool detailed = false) const {
-        if (detailed) {
-            std::string data_str = "[";
-            for (auto &d : data) {
-                data_str += std::to_string(d) + ", ";
-            }
-            data_str += "]";
-            return "x: " + std::to_string(x) + ", y: " + std::to_string(y) +
-                   ", data: " + data_str;
-        }
-        return "x: " + std::to_string(x) + ", y: " + std::to_string(y);
-    }
-};
-
-struct ClusterV2 {
-    ClusterV2_ cluster;
-    int32_t frame_number;
-    std::string to_string() const {
-        return "frame_number: " + std::to_string(frame_number) + ", " +
-               cluster.to_string();
-    }
-};
-
-/**
- * @brief
- * important not: fp always points to the clusters header and does not point to
- * individual clusters
- *
- */
-class ClusterFileV2 {
-    std::filesystem::path m_fpath;
-    std::string m_mode;
-    FILE *fp{nullptr};
-
-    void check_open() {
-        if (!fp)
-            throw std::runtime_error(
-                fmt::format("File: {} not open", m_fpath.string()));
-    }
-
-  public:
-    ClusterFileV2(std::filesystem::path const &fpath, std::string const &mode)
-        : m_fpath(fpath), m_mode(mode) {
-        if (m_mode != "r" && m_mode != "w")
-            throw std::invalid_argument("mode must be 'r' or 'w'");
-        if (m_mode == "r" && !std::filesystem::exists(m_fpath))
-            throw std::invalid_argument("File does not exist");
-        if (mode == "r") {
-            fp = fopen(fpath.string().c_str(), "rb");
-        } else if (mode == "w") {
-            if (std::filesystem::exists(fpath)) {
-                fp = fopen(fpath.string().c_str(), "r+b");
-            } else {
-                fp = fopen(fpath.string().c_str(), "wb");
-            }
-        }
-        if (fp == nullptr) {
-            throw std::runtime_error("Failed to open file");
-        }
-    }
-    ~ClusterFileV2() { close(); }
-    std::vector<ClusterV2> read() {
-        check_open();
-
-        ClusterHeader header;
-        fread(&header, sizeof(ClusterHeader), 1, fp);
-        std::vector<ClusterV2_> clusters_(header.n_clusters);
-        fread(clusters_.data(), sizeof(ClusterV2_), header.n_clusters, fp);
-        std::vector<ClusterV2> clusters;
-        for (auto &c : clusters_) {
-            ClusterV2 cluster;
-            cluster.cluster = std::move(c);
-            cluster.frame_number = header.frame_number;
-            clusters.push_back(cluster);
-        }
-
-        return clusters;
-    }
-    std::vector<std::vector<ClusterV2>> read(int n_frames) {
-        std::vector<std::vector<ClusterV2>> clusters;
-        for (int i = 0; i < n_frames; i++) {
-            clusters.push_back(read());
-        }
-        return clusters;
-    }
-
-    size_t write(std::vector<ClusterV2> const &clusters) {
-        check_open();
-        if (m_mode != "w")
-            throw std::runtime_error("File not opened in write mode");
-        if (clusters.empty())
-            return 0;
-
-        ClusterHeader header;
-        header.frame_number = clusters[0].frame_number;
-        header.n_clusters = clusters.size();
-        fwrite(&header, sizeof(ClusterHeader), 1, fp);
-        for (auto &c : clusters) {
-            fwrite(&c.cluster, sizeof(ClusterV2_), 1, fp);
-        }
-        return clusters.size();
-    }
-
-    size_t write(std::vector<std::vector<ClusterV2>> const &clusters) {
-        check_open();
-        if (m_mode != "w")
-            throw std::runtime_error("File not opened in write mode");
-
-        size_t n_clusters = 0;
-        for (auto &c : clusters) {
-            n_clusters += write(c);
-        }
-        return n_clusters;
-    }
-
-    int seek_to_begin() { return fseek(fp, 0, SEEK_SET); }
-    int seek_to_end() { return fseek(fp, 0, SEEK_END); }
-
-    int32_t frame_number() {
-        auto pos = ftell(fp);
-        ClusterHeader header;
-        fread(&header, sizeof(ClusterHeader), 1, fp);
-        fseek(fp, pos, SEEK_SET);
-        return header.frame_number;
-    }
-
-    void close() {
-        if (fp) {
-            fclose(fp);
-            fp = nullptr;
-        }
-    }
-};
-} // namespace aare

include/aare/ClusterVector.hpp

@@ -18,256 +18,6 @@ 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
- * contiguous memory buffer to store the clusters. It is templated on the data
- * type and the coordinate type of the clusters.
- * @note push_back can invalidate pointers to elements in the container
- * @warning ClusterVector is currently move only to catch unintended copies, but
- * this might change since there are probably use cases where copying is needed.
- * @tparam T data type of the pixels in the cluster
- * @tparam CoordType data type of the x and y coordinates of the cluster
- * (normally int16_t)
- */
-#if 0
-template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType>
-class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
-
-    std::byte *m_data{};
-    size_t m_size{0};
-    size_t m_capacity;
-    uint64_t m_frame_number{0}; // TODO! Check frame number size and type
-    /** 
-    Format string used in the python bindings to create a numpy
-    array from the buffer
-    = - native byte order
-    h - short
-    d - double
-    i - int
-    */
-    constexpr static char m_fmt_base[] = "=h:x:\nh:y:\n({},{}){}:data:";
-
-  public:
-    using value_type = T;
-    using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
-
-    /**
-     * @brief Construct a new ClusterVector object
-     * @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
-     */
-    ClusterVector(size_t capacity = 1024, uint64_t frame_number = 0)
-        : m_capacity(capacity), m_frame_number(frame_number) {
-        allocate_buffer(m_capacity);
-    }
-
-    ~ClusterVector() { delete[] m_data; }
-
-    // Move constructor
-    ClusterVector(ClusterVector &&other) noexcept
-        : m_data(other.m_data), m_size(other.m_size),
-          m_capacity(other.m_capacity), m_frame_number(other.m_frame_number) {
-        other.m_data = nullptr;
-        other.m_size = 0;
-        other.m_capacity = 0;
-    }
-
-    // Move assignment operator
-    ClusterVector &operator=(ClusterVector &&other) noexcept {
-        if (this != &other) {
-            delete[] m_data;
-            m_data = other.m_data;
-            m_size = other.m_size;
-            m_capacity = other.m_capacity;
-            m_frame_number = other.m_frame_number;
-            other.m_data = nullptr;
-            other.m_size = 0;
-            other.m_capacity = 0;
-            other.m_frame_number = 0;
-        }
-        return *this;
-    }
-
-    /**
-     * @brief Reserve space for at least capacity clusters
-     * @param capacity number of clusters to reserve space for
-     * @note If capacity is less than the current capacity, the function does
-     * nothing.
-     */
-    void reserve(size_t capacity) {
-        if (capacity > m_capacity) {
-            allocate_buffer(capacity);
-        }
-    }
-
-    /**
-     * @brief Add a cluster to the vector
-     */
-    void push_back(const ClusterType &cluster) {
-        if (m_size == m_capacity) {
-            allocate_buffer(m_capacity * 2);
-        }
-        std::byte *ptr = element_ptr(m_size);
-        *reinterpret_cast<CoordType *>(ptr) = cluster.x;
-        ptr += sizeof(CoordType);
-        *reinterpret_cast<CoordType *>(ptr) = cluster.y;
-        ptr += sizeof(CoordType);
-
-        std::memcpy(ptr, cluster.data, ClusterSizeX * ClusterSizeY * sizeof(T));
-
-        m_size++;
-    }
-
-    ClusterVector &operator+=(const ClusterVector &other) {
-        if (m_size + other.m_size > m_capacity) {
-            allocate_buffer(m_capacity + other.m_size);
-        }
-        std::copy(other.m_data, other.m_data + other.m_size * item_size(),
-                  m_data + m_size * item_size());
-        m_size += other.m_size;
-        return *this;
-    }
-
-    /**
-     * @brief Sum the pixels in each cluster
-     * @return std::vector<T> vector of sums for each cluster
-     */
-    /*
-    std::vector<T> sum() {
-        std::vector<T> sums(m_size);
-        const size_t stride = item_size();
-        const size_t n_pixels = ClusterSizeX * ClusterSizeY;
-        std::byte *ptr = m_data + 2 * sizeof(CoordType); // skip x and y
-
-        for (size_t i = 0; i < m_size; i++) {
-            sums[i] =
-                std::accumulate(reinterpret_cast<T *>(ptr),
-                                reinterpret_cast<T *>(ptr) + n_pixels, T{});
-            ptr += stride;
-        }
-        return sums;
-    }
-    */
-
-    /**
-     * @brief Sum the pixels in the 2x2 subcluster with the biggest pixel sum in
-     * each cluster
-     * @return std::vector<T> vector of sums for each cluster
-     */ //TODO if underlying container is a vector use std::for_each
-    /*
-    std::vector<T> sum_2x2() {
-        std::vector<T> sums_2x2(m_size);
-
-        for (size_t i = 0; i < m_size; i++) {
-            sums_2x2[i] = at(i).max_sum_2x2;
-        }
-        return sums_2x2;
-    }
-    */
-
-    /**
-     * @brief Return the number of clusters in the vector
-     */
-    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
-     * reallocation.
-     */
-    size_t capacity() const { return m_capacity; }
-
-    /**
-     * @brief Return the size in bytes of a single cluster
-     */
-    size_t item_size() const {
-        return 2 * sizeof(CoordType) + ClusterSizeX * ClusterSizeY * sizeof(T);
-    }
-
-    /**
-     * @brief Return the offset in bytes for the i-th cluster
-     */
-    size_t element_offset(size_t i) const { return item_size() * i; }
-
-    /**
-     * @brief Return a pointer to the i-th cluster
-     */
-    std::byte *element_ptr(size_t i) { return m_data + element_offset(i); }
-
-    /**
-     * @brief Return a pointer to the i-th cluster
-     */
-    const std::byte *element_ptr(size_t i) const {
-        return m_data + element_offset(i);
-    }
-
-    std::byte *data() { return m_data; }
-    std::byte const *data() const { return m_data; }
-
-    /**
-     * @brief Return a reference to the i-th cluster casted to type V
-     * @tparam V type of the cluster
-     */
-    ClusterType &at(size_t i) {
-        return *reinterpret_cast<ClusterType *>(element_ptr(i));
-    }
-
-    const ClusterType &at(size_t i) const {
-        return *reinterpret_cast<const ClusterType *>(element_ptr(i));
-    }
-
-    template <typename V> const V &at(size_t i) const {
-        return *reinterpret_cast<const V *>(element_ptr(i));
-    }
-
-    const std::string_view fmt_base() const {
-        // TODO! how do we match on coord_t?
-        return m_fmt_base;
-    }
-
-    /** 
-     * @brief Return the frame number of the clusters. 0 is used to indicate
-     * that the clusters come from many frames
-     */
-    uint64_t frame_number() const { return m_frame_number; }
-
-    void set_frame_number(uint64_t frame_number) {
-        m_frame_number = frame_number;
-    }
-
-    /** 
-     * @brief Resize the vector to contain new_size clusters. If new_size is
-     * greater than the current capacity, a new buffer is allocated. If the size
-     * is smaller no memory is freed, size is just updated.
-     * @param new_size new size of the vector
-     * @warning The additional clusters are not initialized
-     */
-    void resize(size_t new_size) {
-        // TODO! Should we initialize the new clusters?
-        if (new_size > m_capacity) {
-            allocate_buffer(new_size);
-        }
-        m_size = new_size;
-    }
-
-  private:
-    void allocate_buffer(size_t new_capacity) {
-        size_t num_bytes = item_size() * new_capacity;
-        std::byte *new_data = new std::byte[num_bytes]{};
-        std::copy(m_data, m_data + item_size() * m_size, new_data);
-        delete[] m_data;
-        m_data = new_data;
-        m_capacity = new_capacity;
-    }
-};
-#endif
-
 /**
  * @brief ClusterVector is a container for clusters of various sizes. It
  * uses a contiguous memory buffer to store the clusters. It is templated on
@@ -285,7 +35,7 @@ template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
 class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
 
     std::vector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> m_data{};
-    uint64_t m_frame_number{0}; // TODO! Check frame number size and type
+    int32_t m_frame_number{0}; // TODO! Check frame number size and type
 
   public:
     using value_type = T;
@@ -297,7 +47,7 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
      * @param frame_number frame number of the clusters. Default is 0, which is
      * also used to indicate that the clusters come from many frames
      */
-    ClusterVector(size_t capacity = 300, int32_t frame_number = 0)
+    ClusterVector(size_t capacity = 300, uint64_t frame_number = 0)
         : m_frame_number(frame_number) {
         m_data.reserve(capacity);
     }
@@ -319,6 +69,33 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
         return *this;
     }
 
+    /**
+     * @brief Sum the pixels in each cluster
+     * @return std::vector<T> vector of sums for each cluster
+     */
+    std::vector<T> sum() {
+        std::vector<T> sums(m_data.size());
+
+        for (size_t i = 0; i < m_data.size(); i++) {
+            sums[i] = at(i).sum();
+        }
+        return sums;
+    }
+
+    /**
+     * @brief Sum the pixels in the 2x2 subcluster with the biggest pixel sum in
+     * each cluster
+     * @return std::vector<T> vector of sums for each cluster
+     */ //TODO if underlying container is a vector use std::for_each
+    std::vector<T> sum_2x2() {
+        std::vector<T> sums_2x2(m_data.size());
+
+        for (size_t i = 0; i < m_data.size(); i++) {
+            sums_2x2[i] = at(i).max_sum_2x2().first;
+        }
+        return sums_2x2;
+    }
+
     /**
      * @brief Reserve space for at least capacity clusters
      * @param capacity number of clusters to reserve space for
@@ -361,7 +138,8 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
      * @brief Return the size in bytes of a single cluster
      */
     size_t item_size() const {
-        return 2 * sizeof(CoordType) + ClusterSizeX * ClusterSizeY * sizeof(T);
+        return sizeof(ClusterType); // 2 * sizeof(CoordType) + ClusterSizeX *
+                                    // ClusterSizeY * sizeof(T);
     }
 
     ClusterType *data() { return m_data.data(); }
@@ -379,9 +157,9 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
      * @brief Return the frame number of the clusters. 0 is used to indicate
      * that the clusters come from many frames
      */
-    uint64_t frame_number() const { return m_frame_number; }
+    int32_t frame_number() const { return m_frame_number; }
 
-    void set_frame_number(uint64_t frame_number) {
+    void set_frame_number(int32_t frame_number) {
         m_frame_number = frame_number;
     }
 };

python/CMakeLists.txt

@@ -28,6 +28,9 @@ target_link_libraries(_aare PRIVATE aare_core aare_compiler_flags)
 set( PYTHON_FILES
     aare/__init__.py
     aare/CtbRawFile.py
+    aare/ClusterFinder.py
+    aare/ClusterVector.py
+
     aare/func.py
     aare/RawFile.py
     aare/transform.py
@@ -35,6 +38,7 @@ set( PYTHON_FILES
     aare/utils.py
 )
 
+
 # Copy the python files to the build directory
 foreach(FILE ${PYTHON_FILES})
     configure_file(${FILE} ${CMAKE_BINARY_DIR}/${FILE}  )

python/aare/ClusterFinder.py

@@ -0,0 +1,14 @@
+
+from ._aare import ClusterFinder_Cluster3x3i
+import numpy as np
+
+def ClusterFinder(image_size, cluster_size, n_sigma=5, dtype = np.int32, capacity = 1024):
+    """
+    Factory function to create a ClusterFinder object. Provides a cleaner syntax for 
+    the templated ClusterFinder in C++.
+    """
+    if dtype == np.int32 and cluster_size == (3,3):
+        return ClusterFinder_Cluster3x3i(image_size, n_sigma = n_sigma, capacity=capacity)
+    else:
+        #TODO! add the other formats
+        raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.")

python/aare/ClusterVector.py

@@ -0,0 +1,11 @@
+
+
+from ._aare import ClusterVector_Cluster3x3i
+import numpy as np
+
+def ClusterVector(cluster_size, dtype = np.int32):
+
+    if dtype == np.int32 and cluster_size == (3,3):
+        return ClusterVector_Cluster3x3i()
+    else:
+        raise ValueError(f"Unsupported dtype: {dtype}. Only np.int32 is supported.")

python/aare/__init__.py

@@ -11,8 +11,12 @@ from ._aare import ROI
 
 # from ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i
 
+from .ClusterFinder import ClusterFinder
+from .ClusterVector import ClusterVector
+
 from ._aare import fit_gaus, fit_pol1
 from ._aare import Interpolator
+from ._aare import calculate_eta2
 from .CtbRawFile import CtbRawFile
 from .RawFile import RawFile
 from .ScanParameters import ScanParameters

python/src/bind_ClusterVector.hpp

@@ -0,0 +1,103 @@
+#include "aare/ClusterCollector.hpp"
+#include "aare/ClusterFileSink.hpp"
+#include "aare/ClusterFinder.hpp"
+#include "aare/ClusterFinderMT.hpp"
+#include "aare/ClusterVector.hpp"
+#include "aare/NDView.hpp"
+#include "aare/Pedestal.hpp"
+#include "np_helper.hpp"
+
+#include <cstdint>
+#include <filesystem>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl_bind.h>
+
+namespace py = pybind11;
+using pd_type = double;
+
+using namespace aare;
+
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wunused-parameter"
+
+
+template <typename Type, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = uint16_t>
+void define_ClusterVector(py::module &m, const std::string &typestr) {
+    using ClusterType =
+        Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType, void>;
+    auto class_name = fmt::format("ClusterVector_{}", typestr);
+
+    py::class_<ClusterVector<
+        Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType, void>, void>>(
+        m, class_name.c_str(),
+        py::buffer_protocol())
+
+        .def(py::init()) // TODO change!!!
+
+        .def("push_back",
+             [](ClusterVector<ClusterType> &self, const ClusterType &cluster) {
+                 self.push_back(cluster);
+             })
+
+        .def("sum", [](ClusterVector<ClusterType> &self) {
+            auto *vec = new std::vector<Type>(self.sum());
+            return return_vector(vec);
+        })
+        .def_property_readonly("size", &ClusterVector<ClusterType>::size)
+        .def("item_size", &ClusterVector<ClusterType>::item_size)
+        .def_property_readonly("fmt",
+                               [typestr](ClusterVector<ClusterType> &self) {
+                                   return fmt_format<ClusterType>;
+                               })
+
+        .def_property_readonly("cluster_size_x",
+                               &ClusterVector<ClusterType>::cluster_size_x)
+        .def_property_readonly("cluster_size_y",
+                               &ClusterVector<ClusterType>::cluster_size_y)
+        .def_property_readonly("capacity",
+                               &ClusterVector<ClusterType>::capacity)
+        .def_property("frame_number", &ClusterVector<ClusterType>::frame_number,
+                      &ClusterVector<ClusterType>::set_frame_number)
+        .def_buffer(
+            [typestr](ClusterVector<ClusterType> &self) -> py::buffer_info {
+                return py::buffer_info(
+                    self.data(),             /* Pointer to buffer */
+                    self.item_size(),        /* Size of one scalar */
+                    fmt_format<ClusterType>, /* Format descriptor */
+                    1,                       /* Number of dimensions */
+                    {self.size()},           /* Buffer dimensions */
+                    {self.item_size()} /* Strides (in bytes) for each index */
+                );
+            });
+
+        // Free functions using ClusterVector
+        m.def("hitmap",
+            [](std::array<size_t, 2> image_size, ClusterVector<ClusterType> &cv) {
+                
+                // Create a numpy array to hold the hitmap
+                // The shape of the array is (image_size[0], image_size[1])
+                // note that the python array is passed as [row, col] which
+                // is the opposite of the clusters [x,y]
+                py::array_t<int32_t> hitmap(image_size);
+                auto r = hitmap.mutable_unchecked<2>();
+    
+                // Initialize hitmap to 0
+                for (py::ssize_t i = 0; i < r.shape(0); i++)
+                    for (py::ssize_t j = 0; j < r.shape(1); j++)
+                        r(i, j) = 0;
+    
+
+                // Loop over the clusters and increment the hitmap
+                // Skip out of bound clusters
+                for (const auto& cluster : cv) {
+                    auto x = cluster.x;
+                    auto y = cluster.y;
+                    if(x<image_size[1] && y<image_size[0])
+                        r(cluster.y, cluster.x) += 1;
+                }
+
+                return hitmap;
+            });
+}

python/src/cluster.hpp

@@ -64,53 +64,8 @@ void define_cluster(py::module &m, const std::string &typestr) {
     */
 }
 
-template <typename Type, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
-          typename CoordType = uint16_t>
-void define_cluster_vector(py::module &m, const std::string &typestr) {
-    using ClusterType =
-        Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType, void>;
-    auto class_name = fmt::format("ClusterVector_{}", typestr);
 
-    py::class_<ClusterVector<
-        Cluster<Type, ClusterSizeX, ClusterSizeY, CoordType, void>, void>>(
-        m, class_name.c_str(),
-        py::buffer_protocol())
 
-        .def(py::init()) // TODO change!!!
-
-        .def("push_back",
-             [](ClusterVector<ClusterType> &self, const ClusterType &cluster) {
-                 self.push_back(cluster);
-             })
-
-        // implement push_back
-        .def_property_readonly("size", &ClusterVector<ClusterType>::size)
-        .def("item_size", &ClusterVector<ClusterType>::item_size)
-        .def_property_readonly("fmt",
-                               [typestr](ClusterVector<ClusterType> &self) {
-                                   return fmt_format<ClusterType>;
-                               })
-
-        .def_property_readonly("cluster_size_x",
-                               &ClusterVector<ClusterType>::cluster_size_x)
-        .def_property_readonly("cluster_size_y",
-                               &ClusterVector<ClusterType>::cluster_size_y)
-        .def_property_readonly("capacity",
-                               &ClusterVector<ClusterType>::capacity)
-        .def_property("frame_number", &ClusterVector<ClusterType>::frame_number,
-                      &ClusterVector<ClusterType>::set_frame_number)
-        .def_buffer(
-            [typestr](ClusterVector<ClusterType> &self) -> py::buffer_info {
-                return py::buffer_info(
-                    self.data(),             /* Pointer to buffer */
-                    self.item_size(),        /* Size of one scalar */
-                    fmt_format<ClusterType>, /* Format descriptor */
-                    1,                       /* Number of dimensions */
-                    {self.size()},           /* Buffer dimensions */
-                    {self.item_size()} /* Strides (in bytes) for each index */
-                );
-            });
-}
 
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
           typename CoordType = uint16_t>
@@ -252,25 +207,5 @@ void define_cluster_finder_bindings(py::module &m, const std::string &typestr) {
             },
             py::arg(), py::arg("frame_number") = 0);
 
-    m.def("hitmap",
-          [](std::array<size_t, 2> image_size, ClusterVector<ClusterType> &cv) {
-              py::array_t<int32_t> hitmap(image_size);
-              auto r = hitmap.mutable_unchecked<2>();
-
-              // Initialize hitmap to 0
-              for (py::ssize_t i = 0; i < r.shape(0); i++)
-                  for (py::ssize_t j = 0; j < r.shape(1); j++)
-                      r(i, j) = 0;
-
-              size_t stride = cv.item_size();
-              auto ptr = cv.data();
-              for (size_t i = 0; i < cv.size(); i++) {
-                  auto x = *reinterpret_cast<int16_t *>(ptr);
-                  auto y = *reinterpret_cast<int16_t *>(ptr + sizeof(int16_t));
-                  r(y, x) += 1;
-                  ptr += stride;
-              }
-              return hitmap;
-          });
 }
 #pragma GCC diagnostic pop

python/src/module.cpp

@@ -1,4 +1,9 @@
 // Files with bindings to the different classes
+
+//New style file naming
+#include "bind_ClusterVector.hpp"
+
+//TODO! migrate the other names
 #include "cluster.hpp"
 #include "cluster_file.hpp"
 #include "ctb_raw_file.hpp"
@@ -39,12 +44,12 @@ PYBIND11_MODULE(_aare, m) {
     define_cluster_file_io_bindings<float, 2, 2, uint16_t>(m, "Cluster2x2f");
     define_cluster_file_io_bindings<double, 2, 2, uint16_t>(m, "Cluster2x2d");
 
-    define_cluster_vector<int, 3, 3, uint16_t>(m, "Cluster3x3i");
-    define_cluster_vector<double, 3, 3, uint16_t>(m, "Cluster3x3d");
-    define_cluster_vector<float, 3, 3, uint16_t>(m, "Cluster3x3f");
-    define_cluster_vector<int, 2, 2, uint16_t>(m, "Cluster2x2i");
-    define_cluster_vector<double, 2, 2, uint16_t>(m, "Cluster2x2d");
-    define_cluster_vector<float, 2, 2, uint16_t>(m, "Cluster2x2f");
+    define_ClusterVector<int, 3, 3, uint16_t>(m, "Cluster3x3i");
+    define_ClusterVector<double, 3, 3, uint16_t>(m, "Cluster3x3d");
+    define_ClusterVector<float, 3, 3, uint16_t>(m, "Cluster3x3f");
+    define_ClusterVector<int, 2, 2, uint16_t>(m, "Cluster2x2i");
+    define_ClusterVector<double, 2, 2, uint16_t>(m, "Cluster2x2d");
+    define_ClusterVector<float, 2, 2, uint16_t>(m, "Cluster2x2f");
 
     define_cluster_finder_bindings<int, 3, 3, uint16_t>(m, "Cluster3x3i");
     define_cluster_finder_bindings<double, 3, 3, uint16_t>(m, "Cluster3x3d");

python/tests/test_Cluster.py

@@ -1,12 +1,12 @@
 import pytest 
 import numpy as np
 
-import aare._aare as aare
+from aare import _aare #import the C++ module
 from conftest import test_data_path
 
 
 def test_cluster_vector_can_be_converted_to_numpy():
-    cv = aare.ClusterVector_Cluster3x3i()
+    cv = _aare.ClusterVector_Cluster3x3i()
     arr = np.array(cv, copy=False)
     assert arr.shape == (0,)  # 4 for x, y, size, energy and 9 for the cluster data
 
@@ -14,24 +14,23 @@ def test_cluster_vector_can_be_converted_to_numpy():
 def test_ClusterVector(): 
     """Test ClusterVector""" 
 
-    clustervector = aare.ClusterVector_Cluster3x3i()
+    clustervector = _aare.ClusterVector_Cluster3x3i()
     assert clustervector.cluster_size_x == 3
     assert clustervector.cluster_size_y == 3
     assert clustervector.item_size() == 4+9*4
     assert clustervector.frame_number == 0
-    assert clustervector.capacity == 1024
     assert clustervector.size == 0
 
-    cluster = aare.Cluster3x3i(0,0,np.ones(9, dtype=np.int32))
+    cluster = _aare.Cluster3x3i(0,0,np.ones(9, dtype=np.int32))
     
     clustervector.push_back(cluster) 
     assert clustervector.size == 1
 
     with pytest.raises(TypeError):  # Or use the appropriate exception type
-        clustervector.push_back(aare.Cluster2x2i(0,0,np.ones(4, dtype=np.int32)))
+        clustervector.push_back(_aare.Cluster2x2i(0,0,np.ones(4, dtype=np.int32)))
 
     with pytest.raises(TypeError): 
-        clustervector.push_back(aare.Cluster3x3f(0,0,np.ones(9, dtype=np.float32))) 
+        clustervector.push_back(_aare.Cluster3x3f(0,0,np.ones(9, dtype=np.float32))) 
 
 def test_Interpolator(): 
     """Test Interpolator"""
@@ -41,13 +40,13 @@ def test_Interpolator():
     ybins = np.linspace(0, 5, 30, dtype=np.float64)
 
     etacube = np.zeros(shape=[30, 30, 20], dtype=np.float64)
-    interpolator = aare.Interpolator(etacube, xbins, ybins, ebins)
+    interpolator = _aare.Interpolator(etacube, xbins, ybins, ebins)
 
     assert interpolator.get_ietax().shape == (30,30,20)
     assert interpolator.get_ietay().shape == (30,30,20)
-    clustervector = aare.ClusterVector_Cluster3x3i()
+    clustervector = _aare.ClusterVector_Cluster3x3i()
     
-    cluster = aare.Cluster3x3i(0,0, np.ones(9, dtype=np.int32))
+    cluster = _aare.Cluster3x3i(0,0, np.ones(9, dtype=np.int32))
     clustervector.push_back(cluster) 
 
     interpolated_photons = interpolator.interpolate(clustervector)
@@ -58,9 +57,9 @@ def test_Interpolator():
     assert interpolated_photons[0]["y"] == -1
     assert interpolated_photons[0]["energy"] == 4 #eta_sum = 4, dx, dy = -1,-1 m_ietax = 0, m_ietay = 0
 
-    clustervector = aare.ClusterVector_Cluster2x2i()
+    clustervector = _aare.ClusterVector_Cluster2x2i()
 
-    cluster = aare.Cluster2x2i(0,0, np.ones(4, dtype=np.int32))
+    cluster = _aare.Cluster2x2i(0,0, np.ones(4, dtype=np.int32))
     clustervector.push_back(cluster) 
 
     interpolated_photons = interpolator.interpolate(clustervector)
@@ -71,28 +70,15 @@ def test_Interpolator():
     assert interpolated_photons[0]["y"] == 0
     assert interpolated_photons[0]["energy"] == 4
 
-@pytest.mark.files
-def test_cluster_file(test_data_path): 
-    """Test ClusterFile""" 
-    cluster_file = aare.ClusterFile_Cluster3x3i(test_data_path / "clust/single_frame_97_clustrers.clust") 
-    clustervector = cluster_file.read_clusters(10) #conversion does not work
 
-    cluster_file.close()
-
-    assert clustervector.size == 10
-
-    ###reading with wrong file
-    with pytest.raises(TypeError): 
-        cluster_file = aare.ClusterFile_Cluster2x2i(test_data_path / "clust/single_frame_97_clustrers.clust") 
-        cluster_file.close()
 
 def test_calculate_eta(): 
     """Calculate Eta""" 
-    clusters = aare.ClusterVector_Cluster3x3i() 
-    clusters.push_back(aare.Cluster3x3i(0,0, np.ones(9, dtype=np.int32)))
-    clusters.push_back(aare.Cluster3x3i(0,0, np.array([1,1,1,2,2,2,3,3,3])))
+    clusters = _aare.ClusterVector_Cluster3x3i() 
+    clusters.push_back(_aare.Cluster3x3i(0,0, np.ones(9, dtype=np.int32)))
+    clusters.push_back(_aare.Cluster3x3i(0,0, np.array([1,1,1,2,2,2,3,3,3])))
 
-    eta2 = aare.calculate_eta2(clusters)
+    eta2 = _aare.calculate_eta2(clusters)
 
     assert eta2.shape == (2,2)
     assert eta2[0,0] == 0.5
@@ -103,7 +89,7 @@ def test_calculate_eta():
 def test_cluster_finder(): 
     """Test ClusterFinder""" 
 
-    clusterfinder = aare.ClusterFinder_Cluster3x3i([100,100])
+    clusterfinder = _aare.ClusterFinder_Cluster3x3i([100,100])
 
     #frame = np.random.rand(100,100)
     frame = np.zeros(shape=[100,100])
@@ -115,18 +101,7 @@ def test_cluster_finder():
     assert clusters.size == 0
 
 
-#TODO dont understand behavior
-def test_cluster_collector(): 
-    """Test ClusterCollector"""
-
-    clusterfinder = aare.ClusterFinderMT_Cluster3x3i([100,100]) #TODO: no idea what the data is in InputQueue not zero 
-
-    clustercollector = aare.ClusterCollector_Cluster3x3i(clusterfinder)
-
-    cluster_vectors = clustercollector.steal_clusters()
-
-    assert len(cluster_vectors) == 1 #single thread execution
-    assert cluster_vectors[0].size == 0 #
+
 
 
 

python/tests/test_ClusterFile.py

@@ -0,0 +1,64 @@
+
+import pytest 
+import numpy as np
+import boost_histogram as bh
+import time
+from pathlib import Path
+import pickle
+
+from aare import ClusterFile
+from conftest import test_data_path
+
+@pytest.mark.files
+def test_cluster_file(test_data_path): 
+    """Test ClusterFile""" 
+    f =  ClusterFile(test_data_path / "clust/single_frame_97_clustrers.clust") 
+    cv = f.read_clusters(10) #conversion does not work
+
+
+    assert cv.frame_number == 135
+    assert cv.size == 10
+
+    #Known data 
+    #frame_number, num_clusters   [135] 97
+    #[  1 200] [0 1 2 3 4 5 6 7 8]
+    #[  2 201] [ 9 10 11 12 13 14 15 16 17]
+    #[  3 202] [18 19 20 21 22 23 24 25 26]
+    #[  4 203] [27 28 29 30 31 32 33 34 35]
+    #[  5 204] [36 37 38 39 40 41 42 43 44]
+    #[  6 205] [45 46 47 48 49 50 51 52 53]
+    #[  7 206] [54 55 56 57 58 59 60 61 62]
+    #[  8 207] [63 64 65 66 67 68 69 70 71]
+    #[  9 208] [72 73 74 75 76 77 78 79 80]
+    #[ 10 209] [81 82 83 84 85 86 87 88 89]
+
+    #conversion to numpy array
+    arr = np.array(cv, copy = False)
+    
+    assert arr.size == 10
+    for i in range(10):
+        assert arr[i]['x'] == i+1
+
+@pytest.mark.files
+def test_read_clusters_and_fill_histogram(test_data_path): 
+    # Create the histogram
+    n_bins = 100
+    xmin = -100
+    xmax = 1e4
+    hist_aare = bh.Histogram(bh.axis.Regular(n_bins, xmin, xmax))
+
+    fname = test_data_path / "clust/beam_En700eV_-40deg_300V_10us_d0_f0_100.clust"
+
+    #Read clusters and fill the histogram with pixel values
+    with ClusterFile(fname, chunk_size = 10000) as f:
+        for clusters in f:
+            arr = np.array(clusters, copy = False)
+            hist_aare.fill(arr['data'].flat)
+
+
+    #Load the histogram from the pickle file
+    with open(fname.with_suffix('.pkl'), 'rb') as f:
+        hist_py = pickle.load(f)
+        
+    #Compare the two histograms
+    assert hist_aare == hist_py

python/tests/test_ClusterVector.py

@@ -0,0 +1,54 @@
+import pytest 
+import numpy as np
+import boost_histogram as bh
+import time
+from pathlib import Path
+import pickle
+
+from aare import ClusterFile
+from aare import _aare
+from conftest import test_data_path
+
+
+def test_create_cluster_vector():
+    cv = _aare.ClusterVector_Cluster3x3i()
+    assert cv.cluster_size_x == 3
+    assert cv.cluster_size_y == 3
+    assert cv.size == 0
+
+
+def test_push_back_on_cluster_vector():
+    cv = _aare.ClusterVector_Cluster2x2i()
+    assert cv.cluster_size_x == 2
+    assert cv.cluster_size_y == 2
+    assert cv.size == 0
+
+    cluster = _aare.Cluster2x2i(19, 22, np.ones(4, dtype=np.int32))
+    cv.push_back(cluster)
+    assert cv.size == 1
+
+    arr = np.array(cv, copy=False)
+    assert arr[0]['x'] == 19
+    assert arr[0]['y'] == 22
+
+
+def test_make_a_hitmap_from_cluster_vector():
+    cv = _aare.ClusterVector_Cluster3x3i()
+
+    # Push back 4 clusters with different positions
+    cv.push_back(_aare.Cluster3x3i(0, 0, np.ones(9, dtype=np.int32)))
+    cv.push_back(_aare.Cluster3x3i(1, 1, np.ones(9, dtype=np.int32)))
+    cv.push_back(_aare.Cluster3x3i(1, 1, np.ones(9, dtype=np.int32)))
+    cv.push_back(_aare.Cluster3x3i(2, 2, np.ones(9, dtype=np.int32)))
+
+    ref = np.zeros((5, 5), dtype=np.int32)
+    ref[0,0] = 1
+    ref[1,1] = 2
+    ref[2,2] = 1
+
+
+    img = _aare.hitmap((5,5), cv)
+    # print(img)
+    # print(ref)
+    assert (img == ref).all()
+    

src/CalculateEta.test.cpp

@@ -37,7 +37,17 @@ auto get_test_parameters() {
             Eta2<int>{3. / 5, 4. / 6, 1, 11}));
 }
 
-TEST_CASE("calculate_eta2", "[.eta_calculation]") {
+TEST_CASE("compute_largest_2x2_subcluster", "[eta_calculation]") {
+    auto [cluster, expected_eta] = get_test_parameters();
+
+    auto [sum, index] = std::visit(
+        [](const auto &clustertype) { return clustertype.max_sum_2x2(); },
+        cluster);
+    CHECK(expected_eta.c == index);
+    CHECK(expected_eta.sum == sum);
+}
+
+TEST_CASE("calculate_eta2", "[eta_calculation]") {
 
     auto [cluster, expected_eta] = get_test_parameters();
 
@@ -50,3 +60,69 @@ TEST_CASE("calculate_eta2", "[.eta_calculation]") {
     CHECK(eta.c == expected_eta.c);
     CHECK(eta.sum == expected_eta.sum);
 }
+
+
+//3x3 cluster layout (rotated to match the cBottomLeft enum):
+// 6, 7, 8
+// 3, 4, 5
+// 0, 1, 2
+
+
+TEST_CASE("Calculate eta2 for a 3x3 int32 cluster with the largest 2x2 sum in the bottom left",
+          "[eta_calculation]") {
+
+    // Create a 3x3 cluster
+    Cluster<int32_t, 3, 3> cl;
+    cl.x = 0;
+    cl.y = 0;
+    cl.data[0] = 30;
+    cl.data[1] = 23;
+    cl.data[2] = 5;
+    cl.data[3] = 20;
+    cl.data[4] = 50;
+    cl.data[5] = 3;
+    cl.data[6] = 8;
+    cl.data[7] = 2;
+    cl.data[8] = 3;
+    
+    // 8, 2, 3
+    // 20, 50, 3
+    // 30, 23, 5
+
+    auto eta = calculate_eta2(cl);
+    CHECK(eta.c == corner::cBottomLeft);
+    CHECK(eta.x == 50.0 / (20 + 50)); // 4/(3+4)
+    CHECK(eta.y == 50.0 / (23 + 50)); // 4/(1+4)
+    CHECK(eta.sum == 30+23+20+50);
+
+}
+
+TEST_CASE("Calculate eta2 for a 3x3 int32 cluster with the largest 2x2 sum in the top left",
+    "[eta_calculation]") {
+
+// Create a 3x3 cluster
+Cluster<int32_t, 3, 3> cl;
+cl.x = 0;
+cl.y = 0;
+cl.data[0] = 8;
+cl.data[1] = 12;
+cl.data[2] = 5;
+cl.data[3] = 77;
+cl.data[4] = 80;
+cl.data[5] = 3;
+cl.data[6] = 82;
+cl.data[7] = 91;
+cl.data[8] = 3;
+
+// 82, 91, 3
+// 77, 80, 3
+// 8, 12, 5
+
+auto eta = calculate_eta2(cl);
+CHECK(eta.c == corner::cTopLeft);
+CHECK(eta.x == 80. / (77 + 80)); // 4/(3+4)
+CHECK(eta.y == 91.0 / (91 + 80)); // 7/(7+4)
+CHECK(eta.sum == 77+80+82+91);
+
+}
+

src/Cluster.test.cpp

@@ -26,49 +26,9 @@ TEST_CASE("Correct Instantiation of Cluster and ClusterVector",
     CHECK(not is_cluster_v<int>);
     CHECK(is_cluster_v<Cluster<int, 3, 3>>);
 }
-<<<<<<< Updated upstream
-=======
-
-using ClusterTypes =
-    std::variant<Cluster<int, 2, 2>, Cluster<int, 3, 3>, Cluster<int, 5, 5>,
-                 Cluster<int, 4, 2>, Cluster<int, 2, 3>>;
-
-auto get_test_sum_parameters() {
-    return GENERATE(
-        std::make_tuple(ClusterTypes{Cluster<int, 2, 2>{0, 0, {1, 2, 3, 1}}},
-                        std::make_pair(7, 0)),
-        std::make_tuple(
-            ClusterTypes{Cluster<int, 3, 3>{0, 0, {1, 2, 3, 4, 5, 6, 1, 2, 7}}},
-            std::make_pair(20, 3)),
-        std::make_tuple(ClusterTypes{Cluster<int, 5, 5>{
-                            0, 0, {1, 6, 7, 6, 5, 4, 3, 2, 1, 8, 8, 9, 2,
-                                   1, 4, 5, 6, 7, 8, 4, 1, 1, 1, 1, 1}}},
-                        std::make_pair(28, 8)),
-        std::make_tuple(
-            ClusterTypes{Cluster<int, 4, 2>{0, 0, {1, 4, 7, 2, 5, 6, 4, 3}}},
-            std::make_pair(21, 1)),
-        std::make_tuple(
-            ClusterTypes{Cluster<int, 2, 3>{0, 0, {1, 3, 2, 3, 4, 2}}},
-            std::make_pair(11, 1)));
-}
-
-TEST_CASE("compute_largest_2x2_subcluster", "[.cluster]") {
-    auto [cluster, sum_pair] = get_test_sum_parameters();
-
-    auto sum = std::visit(
-        [](const auto &clustertype) { return clustertype.max_sum_2x2(); },
-        cluster);
-    CHECK(sum_pair.first == sum.first);
-    CHECK(sum_pair.second == sum.second);
-}
 
 TEST_CASE("Test sum of Cluster", "[.cluster]") {
     Cluster<int, 2, 2> cluster{0, 0, {1, 2, 3, 4}};
 
     CHECK(cluster.sum() == 10);
-
-    Cluster<int, 2, 3> cluster2x3{0, 0, {1, 3, 2, 3, 4, 2}};
-
-    CHECK(cluster2x3.sum() == 15);
-}
->>>>>>> Stashed changes
+}

src/ClusterFile.test.cpp

@@ -2,21 +2,29 @@
 #include "test_config.hpp"
 
 #include "aare/defs.hpp"
+#include <algorithm>
 #include <catch2/catch_test_macros.hpp>
 #include <filesystem>
 
 using aare::Cluster;
 using aare::ClusterFile;
+using aare::ClusterVector;
 
-TEST_CASE("Read one frame from a a cluster file", "[.files]") {
+TEST_CASE("Read one frame from a cluster file", "[.files]") {
     // We know that the frame has 97 clusters
     auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust";
     REQUIRE(std::filesystem::exists(fpath));
 
     ClusterFile<Cluster<int32_t, 3, 3>> f(fpath);
     auto clusters = f.read_frame();
-    REQUIRE(clusters.size() == 97);
-    REQUIRE(clusters.frame_number() == 135);
+    CHECK(clusters.size() == 97);
+    CHECK(clusters.frame_number() == 135);
+    CHECK(clusters.at(0).x == 1);
+    CHECK(clusters.at(0).y == 200);
+    int32_t expected_cluster_data[] = {0, 1, 2, 3, 4, 5, 6, 7, 8};
+    CHECK(std::equal(std::begin(clusters.at(0).data),
+                     std::end(clusters.at(0).data),
+                     std::begin(expected_cluster_data)));
 }
 
 TEST_CASE("Read one frame using ROI", "[.files]") {
@@ -43,6 +51,13 @@ TEST_CASE("Read one frame using ROI", "[.files]") {
         REQUIRE(c.y >= roi.ymin);
         REQUIRE(c.y <= roi.ymax);
     }
+
+    CHECK(clusters.at(0).x == 1);
+    CHECK(clusters.at(0).y == 200);
+    int32_t expected_cluster_data[] = {0, 1, 2, 3, 4, 5, 6, 7, 8};
+    CHECK(std::equal(std::begin(clusters.at(0).data),
+                     std::end(clusters.at(0).data),
+                     std::begin(expected_cluster_data)));
 }
 
 TEST_CASE("Read clusters from single frame file", "[.files]") {
@@ -154,6 +169,12 @@ TEST_CASE("Read clusters from single frame file", "[.files]") {
         auto clusters = f.read_clusters(50);
         REQUIRE(clusters.size() == 50);
         REQUIRE(clusters.frame_number() == 135);
+        int32_t expected_cluster_data[] = {0, 1, 2, 3, 4, 5, 6, 7, 8};
+        REQUIRE(clusters.at(0).x == 1);
+        REQUIRE(clusters.at(0).y == 200);
+        CHECK(std::equal(std::begin(clusters.at(0).data),
+                         std::end(clusters.at(0).data),
+                         std::begin(expected_cluster_data)));
     }
     SECTION("Read more clusters than available") {
         ClusterFile<Cluster<int32_t, 3, 3>> f(fpath);
@@ -161,24 +182,169 @@ TEST_CASE("Read clusters from single frame file", "[.files]") {
         auto clusters = f.read_clusters(100);
         REQUIRE(clusters.size() == 97);
         REQUIRE(clusters.frame_number() == 135);
+        int32_t expected_cluster_data[] = {0, 1, 2, 3, 4, 5, 6, 7, 8};
+        REQUIRE(clusters.at(0).x == 1);
+        REQUIRE(clusters.at(0).y == 200);
+        CHECK(std::equal(std::begin(clusters.at(0).data),
+                         std::end(clusters.at(0).data),
+                         std::begin(expected_cluster_data)));
     }
     SECTION("Read all clusters") {
         ClusterFile<Cluster<int32_t, 3, 3>> f(fpath);
         auto clusters = f.read_clusters(97);
         REQUIRE(clusters.size() == 97);
         REQUIRE(clusters.frame_number() == 135);
-
         REQUIRE(clusters.at(0).x == 1);
         REQUIRE(clusters.at(0).y == 200);
+        int32_t expected_cluster_data[] = {0, 1, 2, 3, 4, 5, 6, 7, 8};
+        CHECK(std::equal(std::begin(clusters.at(0).data),
+                         std::end(clusters.at(0).data),
+                         std::begin(expected_cluster_data)));
     }
 }
 
-TEST_CASE("Read clusters", "[.files]") {
-    // beam_En700eV_-40deg_300V_10us_d0_f0_100.clust
-    auto fpath = test_data_path() / "clust" /
-                 "beam_En700eV_-40deg_300V_10us_d0_f0_100.clust";
+TEST_CASE("Read clusters from single frame file with ROI", "[.files]") {
+    auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust";
     REQUIRE(std::filesystem::exists(fpath));
 
     ClusterFile<Cluster<int32_t, 3, 3>> f(fpath);
-    auto clusters = f.read_clusters(500);
+
+    aare::ROI roi;
+    roi.xmin = 0;
+    roi.xmax = 50;
+    roi.ymin = 200;
+    roi.ymax = 249;
+    f.set_roi(roi);
+
+    auto clusters = f.read_clusters(10);
+
+    CHECK(clusters.size() == 10);
+    CHECK(clusters.frame_number() == 135);
+    CHECK(clusters.at(0).x == 1);
+    CHECK(clusters.at(0).y == 200);
+    int32_t expected_cluster_data[] = {0, 1, 2, 3, 4, 5, 6, 7, 8};
+    CHECK(std::equal(std::begin(clusters.at(0).data),
+                     std::end(clusters.at(0).data),
+                     std::begin(expected_cluster_data)));
+}
+
+TEST_CASE("Read cluster from multiple frame file", "[.files]") {
+
+    using ClusterType = Cluster<double, 2, 2>;
+
+    auto fpath =
+        test_data_path() / "clust" / "Two_frames_2x2double_test_clusters.clust";
+
+    REQUIRE(std::filesystem::exists(fpath));
+
+    // Two_frames_2x2double_test_clusters.clust
+    //  frame number, num_clusters 0, 4
+    //[10, 20], {0. ,0., 0., 0.}
+    //[11, 30], {1., 1., 1., 1.}
+    //[12, 40], {2., 2., 2., 2.}
+    //[13, 50], {3., 3., 3., 3.}
+    //  1,4
+    //[10, 20], {4., 4., 4., 4.}
+    //[11, 30], {5., 5., 5., 5.}
+    //[12, 40], {6., 6., 6., 6.}
+    //[13, 50], {7., 7., 7., 7.}
+
+    SECTION("Read clusters from both frames") {
+        ClusterFile<ClusterType> f(fpath);
+        auto clusters = f.read_clusters(2);
+        REQUIRE(clusters.size() == 2);
+        REQUIRE(clusters.frame_number() == 0);
+
+        auto clusters1 = f.read_clusters(3);
+
+        REQUIRE(clusters1.size() == 3);
+        REQUIRE(clusters1.frame_number() == 1);
+    }
+
+    SECTION("Read all clusters") {
+        ClusterFile<ClusterType> f(fpath);
+        auto clusters = f.read_clusters(8);
+        REQUIRE(clusters.size() == 8);
+        REQUIRE(clusters.frame_number() == 1);
+    }
+
+    SECTION("Read clusters from one frame") {
+        ClusterFile<ClusterType> f(fpath);
+        auto clusters = f.read_clusters(2);
+        REQUIRE(clusters.size() == 2);
+        REQUIRE(clusters.frame_number() == 0);
+
+        auto clusters1 = f.read_clusters(1);
+
+        REQUIRE(clusters1.size() == 1);
+        REQUIRE(clusters1.frame_number() == 0);
+    }
+}
+
+TEST_CASE("Write cluster with potential padding", "[.files][.ClusterFile]") {
+
+    using ClusterType = Cluster<double, 3, 3>;
+
+    REQUIRE(std::filesystem::exists(test_data_path() / "clust"));
+
+    auto fpath = test_data_path() / "clust" / "single_frame_2_clusters.clust";
+
+    ClusterFile<ClusterType> file(fpath, 1000, "w");
+
+    ClusterVector<ClusterType> clustervec(2);
+    int16_t coordinate = 5;
+    clustervec.push_back(ClusterType{
+        coordinate, coordinate, {0., 0., 0., 0., 0., 0., 0., 0., 0.}});
+    clustervec.push_back(ClusterType{
+        coordinate, coordinate, {0., 0., 0., 0., 0., 0., 0., 0., 0.}});
+
+    file.write_frame(clustervec);
+
+    file.close();
+
+    file.open("r");
+
+    auto read_cluster_vector = file.read_frame();
+
+    CHECK(read_cluster_vector.size() == 2);
+    CHECK(read_cluster_vector.frame_number() == 0);
+
+    CHECK(read_cluster_vector.at(0).x == clustervec.at(0).x);
+    CHECK(read_cluster_vector.at(0).y == clustervec.at(0).y);
+    CHECK(std::equal(clustervec.at(0).data, clustervec.at(0).data + 9,
+                     read_cluster_vector.at(0).data, [](double a, double b) {
+                         return std::abs(a - b) <
+                                std::numeric_limits<double>::epsilon();
+                     }));
+
+    CHECK(read_cluster_vector.at(1).x == clustervec.at(1).x);
+    CHECK(read_cluster_vector.at(1).y == clustervec.at(1).y);
+    CHECK(std::equal(clustervec.at(1).data, std::end(clustervec.at(1).data),
+                     read_cluster_vector.at(1).data, [](double a, double b) {
+                         return std::abs(a - b) <
+                                std::numeric_limits<double>::epsilon();
+                     }));
+}
+
+TEST_CASE("Read frame and modify cluster data", "[.files][.ClusterFile]") {
+    auto fpath = test_data_path() / "clust" / "single_frame_97_clustrers.clust";
+    REQUIRE(std::filesystem::exists(fpath));
+
+    ClusterFile<Cluster<int32_t, 3, 3>> f(fpath);
+
+    auto clusters = f.read_frame();
+    CHECK(clusters.size() == 97);
+    CHECK(clusters.frame_number() == 135);
+
+    int32_t expected_cluster_data[] = {0, 1, 2, 3, 4, 5, 6, 7, 8};
+    clusters.push_back(
+        Cluster<int32_t, 3, 3>{0, 0, {0, 1, 2, 3, 4, 5, 6, 7, 8}});
+
+    CHECK(clusters.size() == 98);
+    CHECK(clusters.at(0).x == 1);
+    CHECK(clusters.at(0).y == 200);
+
+    CHECK(std::equal(std::begin(clusters.at(0).data),
+                     std::end(clusters.at(0).data),
+                     std::begin(expected_cluster_data)));
 }

src/ClusterVector.test.cpp

@@ -8,15 +8,14 @@
 using aare::Cluster;
 using aare::ClusterVector;
 
-
-TEST_CASE("item_size return the size of the cluster stored"){
+TEST_CASE("item_size return the size of the cluster stored") {
     using C1 = Cluster<int32_t, 2, 2>;
     ClusterVector<C1> cv(4);
     CHECK(cv.item_size() == sizeof(C1));
 
-    //Sanity check 
-    //2*2*4 = 16 bytes of data for the cluster
-    // 2*2 = 4 bytes for the x and y coordinates
+    // Sanity check
+    // 2*2*4 = 16 bytes of data for the cluster
+    //  2*2 = 4 bytes for the x and y coordinates
     REQUIRE(cv.item_size() == 20);
 
     using C2 = Cluster<int32_t, 3, 3>;
@@ -30,8 +29,6 @@ TEST_CASE("item_size return the size of the cluster stored"){
     using C4 = Cluster<char, 10, 5>;
     ClusterVector<C4> cv4(4);
     CHECK(cv4.item_size() == sizeof(C4));
-<<<<<<< Updated upstream
-=======
 
     using C5 = Cluster<int32_t, 2, 3>;
     ClusterVector<C5> cv5(4);
@@ -39,12 +36,11 @@ TEST_CASE("item_size return the size of the cluster stored"){
 
     using C6 = Cluster<double, 5, 5>;
     ClusterVector<C6> cv6(4);
-    CHECK(cv6.item_size() == sizeof(C6)); // double uses padding!!!
+    CHECK(cv6.item_size() == sizeof(C6));
 
     using C7 = Cluster<double, 3, 3>;
     ClusterVector<C7> cv7(4);
     CHECK(cv7.item_size() == sizeof(C7));
->>>>>>> Stashed changes
 }
 
 TEST_CASE("ClusterVector 2x2 int32_t capacity 4, push back then read",
@@ -226,32 +222,6 @@ TEST_CASE("Concatenate two cluster vectors where we need to allocate",
     REQUIRE(ptr[3].y == 17);
 }
 
-TEST_CASE("calculate cluster sum", "[.ClusterVector]") {
-    ClusterVector<Cluster<int32_t, 2, 2>> cv1(2);
-    Cluster<int32_t, 2, 2> c1 = {1, 2, {3, 4, 5, 6}};
-    cv1.push_back(c1);
-    Cluster<int32_t, 2, 2> c2 = {6, 7, {8, 9, 10, 11}};
-    cv1.push_back(c2);
-
-    auto sum1 = cv1.sum();
-
-    std::vector<int32_t> expected_sum1{18, 38};
-
-    CHECK(sum1 == expected_sum1);
-
-    ClusterVector<Cluster<int32_t, 3, 3>> cv2(2);
-    Cluster<int32_t, 3, 3> c3 = {1, 2, {3, 4, 5, 6, 1, 7, 8, 1, 1}};
-    cv2.push_back(c3);
-    Cluster<int32_t, 3, 3> c4 = {6, 7, {8, 9, 10, 11, 13, 5, 12, 2, 4}};
-    cv2.push_back(c4);
-
-    auto sum2 = cv2.sum();
-
-    std::vector<int32_t> expected_sum2{36, 74};
-
-    CHECK(sum2 == expected_sum2);
-}
-
 struct ClusterTestData {
     uint8_t ClusterSizeX;
     uint8_t ClusterSizeY;