added cluster finder

include/aare/ClusterFinder.hpp

@@ -0,0 +1,216 @@
+#pragma once
+#include "aare/Dtype.hpp"
+#include "aare/NDArray.hpp"
+#include "aare/NDView.hpp"
+#include "aare/Pedestal.hpp"
+#include <cstddef>
+
+namespace aare {
+
+/** enum to define the event types */
+enum eventType {
+    PEDESTAL,            /** pedestal */
+    NEIGHBOUR,           /** neighbour i.e. below threshold, but in the cluster of a photon */
+    PHOTON,              /** photon i.e. above threshold */
+    PHOTON_MAX,          /** maximum of a cluster satisfying the photon conditions */
+    NEGATIVE_PEDESTAL,   /** negative value, will not be accounted for as pedestal in order to
+                             avoid drift of the pedestal towards negative values */
+    UNDEFINED_EVENT = -1 /** undefined */
+};
+
+class ClusterFinder {
+  public:
+    ClusterFinder(int cluster_sizeX, int cluster_sizeY, double nSigma = 5.0, double threshold = 0.0)
+        : m_cluster_sizeX(cluster_sizeX), m_cluster_sizeY(cluster_sizeY), m_threshold(threshold), m_nSigma(nSigma) {
+
+        c2 = sqrt((cluster_sizeY + 1) / 2 * (cluster_sizeX + 1) / 2);
+        c3 = sqrt(cluster_sizeX * cluster_sizeY);
+    };
+
+    template <typename FRAME_TYPE, typename PEDESTAL_TYPE>
+    std::vector<Cluster> find_clusters_without_threshold(NDView<FRAME_TYPE, 2> frame, Pedestal<PEDESTAL_TYPE> &pedestal,
+                                                         bool late_update = false) {
+        struct pedestal_update {
+            int x;
+            int y;
+            FRAME_TYPE value;
+        };
+        std::vector<pedestal_update> pedestal_updates;
+
+        std::vector<Cluster> clusters;
+        std::vector<std::vector<eventType>> eventMask;
+        for (int i = 0; i < frame.shape(0); i++) {
+            eventMask.push_back(std::vector<eventType>(frame.shape(1)));
+        }
+        long double val;
+        long double max;
+
+        for (int iy = 0; iy < frame.shape(0); iy++) {
+            for (int ix = 0; ix < frame.shape(1); ix++) {
+                // initialize max and total
+                max = std::numeric_limits<FRAME_TYPE>::min();
+                long double total = 0;
+                eventMask[iy][ix] = PEDESTAL;
+
+                for (short ir = -(m_cluster_sizeY / 2); ir < (m_cluster_sizeY / 2) + 1; ir++) {
+                    for (short ic = -(m_cluster_sizeX / 2); ic < (m_cluster_sizeX / 2) + 1; ic++) {
+                        if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) {
+                            val = frame(iy + ir, ix + ic) - pedestal.mean(iy + ir, ix + ic);
+                            total += val;
+                            if (val > max) {
+                                max = val;
+                            }
+                        }
+                    }
+                }
+                auto rms = pedestal.standard_deviation(iy, ix);
+
+                if (frame(iy, ix) - pedestal.mean(iy, ix) < -m_nSigma * rms) {
+                    eventMask[iy][ix] = NEGATIVE_PEDESTAL;
+                    continue;
+                } else if (max > m_nSigma * rms) {
+                    eventMask[iy][ix] = PHOTON;
+
+                } else if (total > c3 * m_nSigma * rms) {
+                    eventMask[iy][ix] = PHOTON;
+                } else{
+                    if (late_update) {
+                        pedestal_updates.push_back({ix, iy, frame(iy, ix)});
+                    } else {
+                        pedestal.push(iy, ix, frame(iy, ix));
+                    }
+                    continue;
+                }
+                if (eventMask[iy][ix] == PHOTON && (frame(iy, ix) - pedestal.mean(iy, ix)) >= max) {
+                    eventMask[iy][ix] = PHOTON_MAX;
+                    Cluster cluster(m_cluster_sizeX, m_cluster_sizeY, Dtype(typeid(PEDESTAL_TYPE)));
+                    cluster.x = ix;
+                    cluster.y = iy;
+                    short i = 0;
+
+                    for (short ir = -(m_cluster_sizeY / 2); ir < (m_cluster_sizeY / 2) + 1; ir++) {
+                        for (short ic = -(m_cluster_sizeX / 2); ic < (m_cluster_sizeX / 2) + 1; ic++) {
+                            if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) {
+                                PEDESTAL_TYPE tmp = static_cast<PEDESTAL_TYPE>(frame(iy + ir, ix + ic)) -
+                                                    pedestal.mean(iy + ir, ix + ic);
+                                cluster.set<PEDESTAL_TYPE>(i, tmp);
+                                i++;
+                            }
+                        }
+                    }
+                    clusters.push_back(cluster);
+                }
+            }
+        }
+        if (late_update) {
+            for (auto &update : pedestal_updates) {
+                pedestal.push(update.y, update.x, update.value);
+            }
+        }
+        return clusters;
+    }
+
+    template <typename FRAME_TYPE, typename PEDESTAL_TYPE>
+    std::vector<Cluster> find_clusters_with_threshold(NDView<FRAME_TYPE, 2> frame, Pedestal<PEDESTAL_TYPE> &pedestal) {
+        assert(m_threshold > 0);
+        std::vector<Cluster> clusters;
+        std::vector<std::vector<eventType>> eventMask;
+        for (int i = 0; i < frame.shape(0); i++) {
+            eventMask.push_back(std::vector<eventType>(frame.shape(1)));
+        }
+        double tthr, tthr1, tthr2;
+
+        NDArray<FRAME_TYPE, 2> rest({frame.shape(0), frame.shape(1)});
+        NDArray<int, 2> nph({frame.shape(0), frame.shape(1)});
+        // convert to n photons
+        // nph = (frame-pedestal.mean()+0.5*m_threshold)/m_threshold; // can be optimized with expression templates?
+        for (int iy = 0; iy < frame.shape(0); iy++) {
+            for (int ix = 0; ix < frame.shape(1); ix++) {
+                auto val = frame(iy, ix) - pedestal.mean(iy, ix);
+                nph(iy, ix) = (val + 0.5 * m_threshold) / m_threshold;
+                nph(iy, ix) = nph(iy, ix) < 0 ? 0 : nph(iy, ix);
+                rest(iy, ix) = val - nph(iy, ix) * m_threshold;
+            }
+        }
+        // iterate over frame pixels
+        for (int iy = 0; iy < frame.shape(0); iy++) {
+            for (int ix = 0; ix < frame.shape(1); ix++) {
+                eventMask[iy][ix] = PEDESTAL;
+                // initialize max and total
+                FRAME_TYPE max = std::numeric_limits<FRAME_TYPE>::min();
+                long double total = 0;
+                if (rest(iy, ix) <= 0.25 * m_threshold) {
+                    pedestal.push(iy, ix, frame(iy, ix));
+                    continue;
+                }
+                eventMask[iy][ix] = NEIGHBOUR;
+                // iterate over cluster pixels around the current pixel (ix,iy)
+                for (short ir = -(m_cluster_sizeY / 2); ir < (m_cluster_sizeY / 2) + 1; ir++) {
+                    for (short ic = -(m_cluster_sizeX / 2); ic < (m_cluster_sizeX / 2) + 1; ic++) {
+                        if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) {
+                            auto val = frame(iy + ir, ix + ic) - pedestal.mean(iy + ir, ix + ic);
+                            total += val;
+                            if (val > max) {
+                                max = val;
+                            }
+                        }
+                    }
+                }
+
+                auto rms = pedestal.standard_deviation(iy, ix);
+                if (m_nSigma == 0) {
+                    tthr = m_threshold;
+                    tthr1 = m_threshold;
+                    tthr2 = m_threshold;
+                } else {
+                    tthr = m_nSigma * rms;
+                    tthr1 = m_nSigma * rms * c3;
+                    tthr2 = m_nSigma * rms * c2;
+
+                    if (m_threshold > 2 * tthr)
+                        tthr = m_threshold - tthr;
+                    if (m_threshold > 2 * tthr1)
+                        tthr1 = tthr - tthr1;
+                    if (m_threshold > 2 * tthr2)
+                        tthr2 = tthr - tthr2;
+                }
+                if (total > tthr1 || max > tthr) {
+                    eventMask[iy][ix] = PHOTON;
+                    nph(iy, ix) += 1;
+                    rest(iy, ix) -= m_threshold;
+                } else {
+                    pedestal.push(iy, ix, frame(iy, ix));
+                    continue;
+                }
+                if (eventMask[iy][ix] == PHOTON && frame(iy, ix) - pedestal.mean(iy, ix) >= max) {
+                    eventMask[iy][ix] = PHOTON_MAX;
+                    Cluster cluster(m_cluster_sizeX, m_cluster_sizeY, Dtype(typeid(FRAME_TYPE)));
+                    cluster.x = ix;
+                    cluster.y = iy;
+                    short i = 0;
+                    for (short ir = -(m_cluster_sizeY / 2); ir < (m_cluster_sizeY / 2) + 1; ir++) {
+                        for (short ic = -(m_cluster_sizeX / 2); ic < (m_cluster_sizeX / 2) + 1; ic++) {
+                            if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) {
+                                auto tmp = frame(iy + ir, ix + ic) - pedestal.mean(iy + ir, ix + ic);
+                                cluster.set<FRAME_TYPE>(i, tmp);
+                                i++;
+                            }
+                        }
+                    }
+                    clusters.push_back(cluster);
+                }
+            }
+        }
+        return clusters;
+    }
+
+  protected:
+    int m_cluster_sizeX;
+    int m_cluster_sizeY;
+    double m_threshold;
+    double m_nSigma;
+    double c2;
+    double c3;
+};
+
+} // namespace aare

include/aare/Pedestal.hpp

@@ -0,0 +1,121 @@
+#pragma once
+#include "aare/Frame.hpp"
+#include "aare/NDArray.hpp"
+#include "aare/NDView.hpp"
+#include <cstddef>
+
+namespace aare {
+
+template <typename SUM_TYPE = double> class Pedestal {
+  public:
+    Pedestal(uint32_t rows, uint32_t cols, uint32_t n_samples = 1000)
+        : m_rows(rows), m_cols(cols), m_freeze(false), m_samples(n_samples),           m_cur_samples(NDArray<uint32_t, 2>({rows, cols}, 0)),m_sum(NDArray<SUM_TYPE, 2>({rows, cols})),
+ m_sum2(NDArray<SUM_TYPE, 2>({rows, cols})) {
+        assert(rows > 0 && cols > 0 && n_samples > 0);
+        m_sum = 0;
+        m_sum2 = 0;
+    }
+    ~Pedestal() = default;
+
+    NDArray<SUM_TYPE, 2> mean() {
+        NDArray<SUM_TYPE, 2> mean_array({m_rows, m_cols});
+        for (uint32_t i = 0; i < m_rows * m_cols; i++) {
+            mean_array(i / m_cols, i % m_cols) = mean(i / m_cols, i % m_cols);
+        }
+        return mean_array;
+    }
+
+    NDArray<SUM_TYPE, 2> variance() {
+        NDArray<SUM_TYPE, 2> variance_array({m_rows, m_cols});
+        for (uint32_t i = 0; i < m_rows * m_cols; i++) {
+            variance_array(i / m_cols, i % m_cols) = variance(i / m_cols, i % m_cols);
+        }
+        return variance_array;
+    }
+
+    NDArray<SUM_TYPE, 2> standard_deviation() {
+        NDArray<SUM_TYPE, 2> standard_deviation_array({m_rows, m_cols});
+        for (uint32_t i = 0; i < m_rows * m_cols; i++) {
+            standard_deviation_array(i / m_cols, i % m_cols) = standard_deviation(i / m_cols, i % m_cols);
+        }
+
+        return standard_deviation_array;
+    }
+    void clear() {
+        for (uint32_t i = 0; i < m_rows * m_cols; i++) {
+            clear(i / m_cols, i % m_cols);
+        }
+    }
+
+    /*
+     * index level operations
+     */
+    SUM_TYPE mean(const uint32_t row, const uint32_t col) const {
+        if (m_cur_samples(row, col) == 0) {
+            return 0.0;
+        }
+        return m_sum(row, col) / m_cur_samples(row, col);
+    }
+    SUM_TYPE variance(const uint32_t row, const uint32_t col) const {
+        if (m_cur_samples(row, col) == 0) {
+            return 0.0;
+        }
+        return m_sum2(row, col) / m_cur_samples(row, col) - mean(row, col) * mean(row, col);
+    }
+    SUM_TYPE standard_deviation(const uint32_t row, const uint32_t col) const { return std::sqrt(variance(row, col)); }
+
+    void clear(const uint32_t row, const uint32_t col) {
+        m_sum(row, col) = 0;
+        m_sum2(row, col) = 0;
+        m_cur_samples(row, col) = 0;
+    }
+    // frame level operations
+    template <typename T> void push(NDView<T, 2> frame) {
+        assert(frame.size() == m_rows * m_cols);
+        // TODO: test the effect of #pragma omp parallel for
+        for (uint32_t index = 0; index < m_rows * m_cols; index++) {
+            push<T>(index / m_cols, index % m_cols, frame(index));
+        }
+    }
+    template <typename T> void push(Frame &frame) {
+        assert(frame.rows() == static_cast<size_t>(m_rows) && frame.cols() == static_cast<size_t>(m_cols));
+        push<T>(frame.view<T>());
+    }
+
+    // getter functions
+    inline uint32_t rows() const { return m_rows; }
+    inline uint32_t cols() const { return m_cols; }
+    inline uint32_t n_samples() const { return m_samples; }
+    inline NDArray<uint32_t, 2> cur_samples() const { return m_cur_samples; }
+    inline NDArray<SUM_TYPE, 2> get_sum() const { return m_sum; }
+    inline NDArray<SUM_TYPE, 2> get_sum2() const { return m_sum2; }
+
+    // pixel level operations (should be refactored to allow users to implement their own pixel level operations)
+    template <typename T> inline void push(const uint32_t row, const uint32_t col, const T val_) {
+        if (m_freeze) {
+            return;
+        }
+        SUM_TYPE val = static_cast<SUM_TYPE>(val_);
+        const uint32_t idx = index(row, col);
+        if (m_cur_samples(idx) < m_samples) {
+            m_sum(idx) += val;
+            m_sum2(idx) += val * val;
+            m_cur_samples(idx)++;
+        } else {
+            m_sum(idx) += val - m_sum(idx) / m_cur_samples(idx);
+            m_sum2(idx) += val * val - m_sum2(idx) / m_cur_samples(idx);
+        }
+    }
+    inline uint32_t index(const uint32_t row, const uint32_t col) const { return row * m_cols + col; };
+    void set_freeze(bool freeze) { m_freeze = freeze; }
+
+  private:
+    uint32_t m_rows;
+    uint32_t m_cols;
+    bool m_freeze;
+    uint32_t m_samples;
+    NDArray<uint32_t, 2> m_cur_samples;
+    NDArray<SUM_TYPE, 2> m_sum;
+    NDArray<SUM_TYPE, 2> m_sum2;
+};
+} // namespace aare

include/aare/VariableSizeClusterFinder.hpp

@@ -0,0 +1,307 @@
+#pragma once
+
+#include <algorithm>
+#include <map>
+#include <unordered_map>
+#include <vector>
+
+#include "aare/NDArray.hpp"
+
+const int MAX_CLUSTER_SIZE = 200;
+namespace aare {
+
+template <typename T> class ClusterFinder {
+  public:
+    struct Hit {
+        int16_t size{};
+        int16_t row{};
+        int16_t col{};
+        uint16_t reserved{}; // for alignment
+        T energy{};
+        T max{};
+
+        // std::vector<int16_t> rows{};
+        // std::vector<int16_t> cols{};
+        int16_t rows[MAX_CLUSTER_SIZE] = {0};
+        int16_t cols[MAX_CLUSTER_SIZE] = {0};
+        double enes[MAX_CLUSTER_SIZE] = {0};
+    };
+
+  private:
+    const std::array<int64_t, 2> shape_;
+    NDView<T, 2> original_;
+    NDArray<int, 2> labeled_;
+    NDArray<int, 2> peripheral_labeled_;
+    NDArray<bool, 2> binary_; // over threshold flag
+    T threshold_;
+    NDView<T, 2> noiseMap;
+    bool use_noise_map = false;
+    int peripheralThresholdFactor_ = 5;
+    int current_label;
+    const std::array<int, 4> di{{0, -1, -1, -1}};              // row ### 8-neighbour by scaning from left to right
+    const std::array<int, 4> dj{{-1, -1, 0, 1}};               // col ### 8-neighbour by scaning from top to bottom
+    const std::array<int, 8> di_{{0, 0, -1, 1, -1, 1, -1, 1}}; // row
+    const std::array<int, 8> dj_{{-1, 1, 0, 0, 1, -1, -1, 1}}; // col
+    std::map<int, int> child;                                  // heirachy: key: child; val: parent
+    std::unordered_map<int, Hit> h_size;
+    std::vector<Hit> hits;
+    // std::vector<std::vector<int16_t>> row
+    int check_neighbours(int i, int j);
+
+  public:
+    ClusterFinder(image_shape shape, T threshold)
+        : shape_(shape), labeled_(shape, 0), peripheral_labeled_(shape, 0), binary_(shape), threshold_(threshold) {
+        hits.reserve(2000);
+    }
+
+    NDArray<int, 2> labeled() { return labeled_; }
+
+    void set_noiseMap(NDView<T, 2> noise_map) {
+        noiseMap = noise_map;
+        use_noise_map = true;
+    }
+    void set_peripheralThresholdFactor(int factor) { peripheralThresholdFactor_ = factor; }
+    void find_clusters(NDView<T, 2> img);
+    void find_clusters_X(NDView<T, 2> img);
+    void rec_FillHit(int clusterIndex, int i, int j);
+    void single_pass(NDView<T, 2> img);
+    void first_pass();
+    void second_pass();
+    void store_clusters();
+
+    std::vector<Hit> steal_hits() {
+        std::vector<Hit> tmp;
+        std::swap(tmp, hits);
+        return tmp;
+    };
+    void clear_hits() { hits.clear(); };
+
+    void print_connections() {
+        fmt::print("Connections:\n");
+        for (auto it = child.begin(); it != child.end(); ++it) {
+            fmt::print("{} -> {}\n", it->first, it->second);
+        }
+    }
+    size_t total_clusters() const {
+        // TODO! fix for stealing
+        return hits.size();
+    }
+
+  private:
+    void add_link(int from, int to) {
+        // we want to add key from -> value to
+        // fmt::print("add_link({},{})\n", from, to);
+        auto it = child.find(from);
+        if (it == child.end()) {
+            child[from] = to;
+        } else {
+            // found need to disambiguate
+            if (it->second == to)
+                return;
+            else {
+                if (it->second > to) {
+                    // child[from] = to;
+                    auto old = it->second;
+                    it->second = to;
+                    add_link(old, to);
+                } else {
+                    // found value is smaller than what we want to link
+                    add_link(to, it->second);
+                }
+            }
+        }
+    }
+};
+template <typename T> int ClusterFinder<T>::check_neighbours(int i, int j) {
+    std::vector<int> neighbour_labels;
+
+    for (int k = 0; k < 4; ++k) {
+        const auto row = i + di[k];
+        const auto col = j + dj[k];
+        if (row >= 0 && col >= 0 && row < shape_[0] && col < shape_[1]) {
+            auto tmp = labeled_.value(i + di[k], j + dj[k]);
+            if (tmp != 0)
+                neighbour_labels.push_back(tmp);
+        }
+    }
+
+    if (neighbour_labels.size() == 0) {
+        return 0;
+    } else {
+
+        // need to sort and add to union field
+        std::sort(neighbour_labels.rbegin(), neighbour_labels.rend());
+        auto first = neighbour_labels.begin();
+        auto last = std::unique(first, neighbour_labels.end());
+        if (last - first == 1)
+            return *neighbour_labels.begin();
+
+        for (auto current = first; current != last - 1; ++current) {
+            auto next = current + 1;
+            add_link(*current, *next);
+        }
+        return neighbour_labels.back(); // already sorted
+    }
+}
+
+template <typename T> void ClusterFinder<T>::find_clusters(NDView<T, 2> img) {
+    original_ = img;
+    labeled_ = 0;
+    peripheral_labeled_ = 0;
+    current_label = 0;
+    child.clear();
+    first_pass();
+    // print_connections();
+    second_pass();
+    store_clusters();
+}
+
+template <typename T> void ClusterFinder<T>::find_clusters_X(NDView<T, 2> img) {
+    original_ = img;
+    int clusterIndex = 0;
+    for (int i = 0; i < shape_[0]; ++i) {
+        for (int j = 0; j < shape_[1]; ++j) {
+            if (use_noise_map)
+                threshold_ = 5 * noiseMap(i, j);
+            if (original_(i, j) > threshold_) {
+                // printf("========== Cluster index: %d\n", clusterIndex);
+                rec_FillHit(clusterIndex, i, j);
+                clusterIndex++;
+            }
+        }
+    }
+    for (const auto &h : h_size)
+        hits.push_back(h.second);
+    h_size.clear();
+}
+
+template <typename T> void ClusterFinder<T>::rec_FillHit(int clusterIndex, int i, int j) {
+    // printf("original_(%d, %d)=%f\n", i, j, original_(i,j));
+    // printf("h_size[%d].size=%d\n", clusterIndex, h_size[clusterIndex].size);
+    if (h_size[clusterIndex].size < MAX_CLUSTER_SIZE) {
+        h_size[clusterIndex].rows[h_size[clusterIndex].size] = i;
+        h_size[clusterIndex].cols[h_size[clusterIndex].size] = j;
+        h_size[clusterIndex].enes[h_size[clusterIndex].size] = original_(i, j);
+    }
+    h_size[clusterIndex].size += 1;
+    h_size[clusterIndex].energy += original_(i, j);
+    if (h_size[clusterIndex].max < original_(i, j)) {
+        h_size[clusterIndex].row = i;
+        h_size[clusterIndex].col = j;
+        h_size[clusterIndex].max = original_(i, j);
+    }
+    original_(i, j) = 0;
+
+    for (int k = 0; k < 8; ++k) { // 8 for 8-neighbour
+        const auto row = i + di_[k];
+        const auto col = j + dj_[k];
+        if (row >= 0 && col >= 0 && row < shape_[0] && col < shape_[1]) {
+            if (use_noise_map)
+                threshold_ = peripheralThresholdFactor_ * noiseMap(row, col);
+            if (original_(row, col) > threshold_) {
+                rec_FillHit(clusterIndex, row, col);
+            } else {
+                // if (h_size[clusterIndex].size < MAX_CLUSTER_SIZE){
+                //     h_size[clusterIndex].size += 1;
+                //     h_size[clusterIndex].rows[h_size[clusterIndex].size] = row;
+                //     h_size[clusterIndex].cols[h_size[clusterIndex].size] = col;
+                //     h_size[clusterIndex].enes[h_size[clusterIndex].size] = original_(row, col);
+                // }// ? weather to include peripheral pixels
+                original_(row, col) = 0; // remove peripheral pixels, to avoid potential influence for pedestal updating
+            }
+        }
+    }
+}
+
+template <typename T> void ClusterFinder<T>::single_pass(NDView<T, 2> img) {
+    original_ = img;
+    labeled_ = 0;
+    current_label = 0;
+    child.clear();
+    first_pass();
+    // print_connections();
+    // second_pass();
+    // store_clusters();
+}
+
+template <typename T> void ClusterFinder<T>::first_pass() {
+
+    for (int i = 0; i < original_.size(); ++i) {
+        if (use_noise_map)
+            threshold_ = 5 * noiseMap(i);
+        binary_(i) = (original_(i) > threshold_);
+    }
+
+    for (int i = 0; i < shape_[0]; ++i) {
+        for (int j = 0; j < shape_[1]; ++j) {
+
+            // do we have someting to process?
+            if (binary_(i, j)) {
+                auto tmp = check_neighbours(i, j);
+                if (tmp != 0) {
+                    labeled_(i, j) = tmp;
+                } else {
+                    labeled_(i, j) = ++current_label;
+                }
+            }
+        }
+    }
+}
+
+template <typename T> void ClusterFinder<T>::second_pass() {
+
+    for (int64_t i = 0; i != labeled_.size(); ++i) {
+        auto current_label = labeled_(i);
+        if (current_label != 0) {
+            auto it = child.find(current_label);
+            while (it != child.end()) {
+                current_label = it->second;
+                it = child.find(current_label);
+                // do this once before doing the second pass?
+                // all values point to the final one...
+            }
+            labeled_(i) = current_label;
+        }
+    }
+}
+
+template <typename T> void ClusterFinder<T>::store_clusters() {
+
+    // Accumulate hit information in a map
+    // Do we always have monotonic increasing
+    // labels? Then vector?
+    // here the translation is label -> Hit
+    std::unordered_map<int, Hit> h_size;
+    for (int i = 0; i < shape_[0]; ++i) {
+        for (int j = 0; j < shape_[1]; ++j) {
+            if (labeled_(i, j) != 0 || false
+                // (i-1 >= 0 and labeled_(i-1, j) != 0) or // another circle of peripheral pixels
+                // (j-1 >= 0 and labeled_(i, j-1) != 0) or
+                // (i+1 < shape_[0] and labeled_(i+1, j) != 0) or
+                // (j+1 < shape_[1] and labeled_(i, j+1) != 0)
+            ) {
+                Hit &record = h_size[labeled_(i, j)];
+                if (record.size < MAX_CLUSTER_SIZE) {
+                    record.rows[record.size] = i;
+                    record.cols[record.size] = j;
+                    record.enes[record.size] = original_(i, j);
+                } else {
+                    continue;
+                }
+                record.size += 1;
+                record.energy += original_(i, j);
+
+                if (record.max < original_(i, j)) {
+                    record.row = i;
+                    record.col = j;
+                    record.max = original_(i, j);
+                }
+            }
+        }
+    }
+
+    for (const auto &h : h_size)
+        hits.push_back(h.second);
+}
+
+} // namespace aare

src/CMakeLists.txt

@@ -23,6 +23,8 @@ if(AARE_TESTS)
             # ${CMAKE_CURRENT_SOURCE_DIR}/test/ProducerConsumerQueue.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/NDArray.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/NDView.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/ClusterFinder.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/Pedestal.test.cpp
             # ${CMAKE_CURRENT_SOURCE_DIR}/test/CircularFifo.test.cpp
             # ${CMAKE_CURRENT_SOURCE_DIR}/test/wrappers.test.cpp
             # ${CMAKE_CURRENT_SOURCE_DIR}/test/Transforms.test.cpp

src/ClusterFinder.test.cpp

@@ -0,0 +1,59 @@
+#include "aare/ClusterFinder.hpp"
+#include "aare/Pedestal.hpp"
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <chrono>
+#include <random>
+
+using namespace aare;
+
+class ClusterFinderUnitTest : public ClusterFinder {
+  public:
+    ClusterFinderUnitTest(int cluster_sizeX, int cluster_sizeY, double nSigma = 5.0, double threshold = 0.0)
+        : ClusterFinder(cluster_sizeX, cluster_sizeY, nSigma, threshold) {}
+    double get_c2() { return c2; }
+    double get_c3() { return c3; }
+    auto get_threshold() { return m_threshold; }
+    auto get_nSigma() { return m_nSigma; }
+    auto get_cluster_sizeX() { return m_cluster_sizeX; }
+    auto get_cluster_sizeY() { return m_cluster_sizeY; }
+};
+
+TEST_CASE("test ClusterFinder constructor") {
+    ClusterFinderUnitTest cf(55, 100);
+    REQUIRE(cf.get_cluster_sizeX() == 55);
+    REQUIRE(cf.get_cluster_sizeY() == 100);
+    REQUIRE(cf.get_threshold() == 0.0);
+    REQUIRE(cf.get_nSigma() == 5.0);
+    double c2 = sqrt((100 + 1) / 2 * (55 + 1) / 2);
+    double c3 = sqrt(55 * 100);
+    // REQUIRE(compare_floats<double>(cf.get_c2(), c2));
+    // REQUIRE(compare_floats<double>(cf.get_c3(), c3));
+    REQUIRE_THAT(cf.get_c2(), Catch::Matchers::WithinRel(c2, 1e-9));
+    REQUIRE_THAT(cf.get_c3(), Catch::Matchers::WithinRel(c3, 1e-9));
+}
+
+TEST_CASE("test cluster finder") {
+    aare::Pedestal pedestal(10, 10, 5);
+    NDArray<double, 2> frame({10, 10});
+    frame = 0;
+    ClusterFinder clusterFinder(3, 3, 1, 1); // 3x3 cluster, 1 nSigma, 1 threshold
+
+    auto clusters = clusterFinder.find_clusters_without_threshold(frame.span(), pedestal);
+
+    REQUIRE(clusters.size() == 0);
+
+    frame(5, 5) = 10;
+    clusters = clusterFinder.find_clusters_without_threshold(frame.span(), pedestal);
+    REQUIRE(clusters.size() == 1);
+    REQUIRE(clusters[0].x == 5);
+    REQUIRE(clusters[0].y == 5);
+    for (int i = 0; i < 3; i++) {
+        for (int j = 0; j < 3; j++) {
+            if (i == 1 && j == 1)
+                REQUIRE(clusters[0].get<double>(i * 3 + j) == 10);
+            else
+                REQUIRE(clusters[0].get<double>(i * 3 + j) == 0);
+        }
+    }
+}

src/Pedestal.test.cpp

@@ -0,0 +1,103 @@
+#include "aare/Pedestal.hpp"
+
+
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <chrono>
+#include <random>
+
+using namespace aare;
+TEST_CASE("test pedestal constructor") {
+    aare::Pedestal pedestal(10, 10, 5);
+    REQUIRE(pedestal.rows() == 10);
+    REQUIRE(pedestal.cols() == 10);
+    REQUIRE(pedestal.n_samples() == 5);
+    for (int i = 0; i < 10; i++) {
+        for (int j = 0; j < 10; j++) {
+            REQUIRE(pedestal.get_sum()(i, j) == 0);
+            REQUIRE(pedestal.get_sum2()(i, j) == 0);
+            REQUIRE(pedestal.cur_samples()(i, j) == 0);
+        }
+    }
+}
+
+TEST_CASE("test pedestal push") {
+    aare::Pedestal pedestal(10, 10, 5);
+    aare::Frame frame(10, 10, Dtype::UINT16);
+    for (int i = 0; i < 10; i++) {
+        for (int j = 0; j < 10; j++) {
+            frame.set<uint16_t>(i, j, i + j);
+        }
+    }
+
+    // test single push
+    pedestal.push<uint16_t>(frame);
+    for (int i = 0; i < 10; i++) {
+        for (int j = 0; j < 10; j++) {
+            REQUIRE(pedestal.get_sum()(i, j) == i + j);
+            REQUIRE(pedestal.get_sum2()(i, j) == (i + j) * (i + j));
+            REQUIRE(pedestal.cur_samples()(i, j) == 1);
+        }
+    }
+
+    // test clear
+    pedestal.clear();
+    for (int i = 0; i < 10; i++) {
+        for (int j = 0; j < 10; j++) {
+            REQUIRE(pedestal.get_sum()(i, j) == 0);
+            REQUIRE(pedestal.get_sum2()(i, j) == 0);
+            REQUIRE(pedestal.cur_samples()(i, j) == 0);
+        }
+    }
+
+    // test number of samples after multiple push
+    for (uint32_t k = 0; k < 50; k++) {
+        pedestal.push<uint16_t>(frame);
+        for (uint32_t i = 0; i < 10; i++) {
+            for (uint32_t j = 0; j < 10; j++) {
+                if (k < 5) {
+                    REQUIRE(pedestal.cur_samples()(i, j) == k + 1);
+                    REQUIRE(pedestal.get_sum()(i, j) == (k + 1) * (i + j));
+                    REQUIRE(pedestal.get_sum2()(i, j) == (k + 1) * (i + j) * (i + j));
+                } else {
+                    REQUIRE(pedestal.cur_samples()(i, j) == 5);
+                    REQUIRE(pedestal.get_sum()(i, j) == 5 * (i + j));
+                    REQUIRE(pedestal.get_sum2()(i, j) == 5 * (i + j) * (i + j));
+                }
+                REQUIRE(pedestal.mean(i, j) == (i + j));
+                REQUIRE(pedestal.variance(i, j) == 0);
+                REQUIRE(pedestal.standard_deviation(i, j) == 0);
+            }
+        }
+    }
+}
+
+TEST_CASE("test pedestal with normal distribution") {
+    const double MEAN = 5.0, STD = 2.0, VAR = STD * STD, TOLERANCE = 0.1;
+
+    unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
+    std::default_random_engine generator(seed);
+    std::normal_distribution<double> distribution(MEAN, STD);
+
+    aare::Pedestal pedestal(3, 5, 10000);
+    for (int i = 0; i < 10000; i++) {
+        aare::Frame frame(3, 5, Dtype::DOUBLE);
+        for (int j = 0; j < 3; j++) {
+            for (int k = 0; k < 5; k++) {
+                frame.set<double>(j, k, distribution(generator));
+            }
+        }
+        pedestal.push<double>(frame);
+    }
+    auto mean = pedestal.mean();
+    auto variance = pedestal.variance();
+    auto standard_deviation = pedestal.standard_deviation();
+
+    for (int i = 0; i < 3; i++) {
+        for (int j = 0; j < 5; j++) {
+            REQUIRE_THAT(mean(i, j), Catch::Matchers::WithinAbs(MEAN, MEAN * TOLERANCE));
+            REQUIRE_THAT(variance(i, j), Catch::Matchers::WithinAbs(VAR, VAR * TOLERANCE));
+            REQUIRE_THAT(standard_deviation(i, j), Catch::Matchers::WithinAbs(STD, STD * TOLERANCE));
+        }
+    }
+}

tests/test_config.hpp.in

@@ -1,7 +1,12 @@
 #pragma once
 #include <filesystem>
+#include <cstdlib>
 
-static constexpr auto test_data_path_str = "@TEST_FILE_PATH@";
-inline auto test_data_path() {
-    return std::filesystem::path(test_data_path_str);
+
+inline auto test_data_path(){
+    if(const char* env_p = std::getenv("AARE_TEST_DATA")){
+        return std::filesystem::path(env_p);
+    }else{
+        throw std::runtime_error("AARE_TEST_DATA_PATH not set");
+    }
 }