WIP

include/aare/ClusterFile.hpp

@@ -1,4 +1,4 @@
-
+#pragma once
 
 #include "aare/defs.hpp"
 #include <filesystem>

include/aare/ClusterFinder.hpp

@@ -1,4 +1,6 @@
 #pragma once
+#include "aare/ClusterFile.hpp"
+#include "aare/ClusterVector.hpp"
 #include "aare/Dtype.hpp"
 #include "aare/NDArray.hpp"
 #include "aare/NDView.hpp"
@@ -9,7 +11,7 @@
 namespace aare {
 
 /** enum to define the event types */
-enum eventType {
+enum class eventType {
     PEDESTAL,   /** pedestal */
     NEIGHBOUR,  /** neighbour i.e. below threshold, but in the cluster of a
                    photon */
@@ -33,118 +35,101 @@ class ClusterFinder {
     Pedestal<PEDESTAL_TYPE> m_pedestal;
 
   public:
-    ClusterFinder(Shape<2> image_size, Shape<2>cluster_size, double nSigma = 5.0,
+    ClusterFinder(Shape<2> image_size, Shape<2> cluster_size,
-                  double threshold = 0.0)
+                  double nSigma = 5.0, double threshold = 0.0)
-        : m_image_size(image_size), m_cluster_sizeX(cluster_size[0]), m_cluster_sizeY(cluster_size[1]),
+        : m_image_size(image_size), m_cluster_sizeX(cluster_size[0]),
-          m_threshold(threshold), m_nSigma(nSigma),
+          m_cluster_sizeY(cluster_size[1]), m_threshold(threshold),
+          m_nSigma(nSigma),
           c2(sqrt((m_cluster_sizeY + 1) / 2 * (m_cluster_sizeX + 1) / 2)),
           c3(sqrt(m_cluster_sizeX * m_cluster_sizeY)),
           m_pedestal(image_size[0], image_size[1]) {
-
+        fmt::print("TypeIndex: {}\n", sizeof(Dtype));
-              // c2 = sqrt((cluster_sizeY + 1) / 2 * (cluster_sizeX + 1) / 2);
+    };
-              // c3 = sqrt(cluster_sizeX * cluster_sizeY);
-          };
 
     void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {
         m_pedestal.push(frame);
     }
 
-    NDArray<PEDESTAL_TYPE, 2> pedestal() {
+    NDArray<PEDESTAL_TYPE, 2> pedestal() { return m_pedestal.mean(); }
-        return m_pedestal.mean();
-    }
 
-    std::vector<DynamicCluster>
+    NDArray<PEDESTAL_TYPE, 2> noise() { return m_pedestal.std(); }
-    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<DynamicCluster> clusters;
+    ClusterVector<PEDESTAL_TYPE>
-        std::vector<std::vector<eventType>> eventMask;
+    find_clusters_without_threshold(NDView<FRAME_TYPE, 2> frame) {
-        for (int i = 0; i < frame.shape(0); i++) {
+        // std::vector<DynamicCluster> clusters;
-            eventMask.push_back(std::vector<eventType>(frame.shape(1)));
+        // std::vector<Cluster> clusters; //Hard coded 3x3 cluster
-        }
+        // clusters.reserve(2000);
-        long double val;
+        ClusterVector<PEDESTAL_TYPE> clusters(m_cluster_sizeX, m_cluster_sizeY);
-        long double max;
+        eventType event_type = eventType::PEDESTAL;
+
+        // TODO! deal with even size clusters
+        // currently 3,3 -> +/- 1
+        //  4,4 -> +/- 2
+        short dy = m_cluster_sizeY / 2;
+        short dx = m_cluster_sizeX / 2;
 
         for (int iy = 0; iy < frame.shape(0); iy++) {
             for (int ix = 0; ix < frame.shape(1); ix++) {
-                // initialize max and total
+                PEDESTAL_TYPE max = std::numeric_limits<FRAME_TYPE>::min();
-                max = std::numeric_limits<FRAME_TYPE>::min();
+                PEDESTAL_TYPE total = 0;
-                long double total = 0;
-                eventMask[iy][ix] = PEDESTAL;
 
-                for (short ir = -(m_cluster_sizeY / 2);
+                for (short ir = -dy; ir < dy + 1; ir++) {
-                     ir < (m_cluster_sizeY / 2) + 1; ir++) {
+                    for (short ic = -dx; ic < dx + 1; ic++) {
-                    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_TYPE val =
-                                  m_pedestal.mean(iy + ir, ix + ic);
+                                frame(iy + ir, ix + ic) -
+                                m_pedestal.mean(iy + ir, ix + ic);
+
                             total += val;
-                            if (val > max) {
+                            max = std::max(max, val);
-                                max = val;
-                            }
                         }
                     }
                 }
-                auto rms = m_pedestal.std(iy, ix);
+                PEDESTAL_TYPE rms = m_pedestal.std(iy, ix);
+                PEDESTAL_TYPE value = (frame(iy, ix) - m_pedestal.mean(iy, ix));
 
-                if (frame(iy, ix) - m_pedestal.mean(iy, ix) < -m_nSigma * rms) {
+                if (value < -m_nSigma * rms) {
-                    eventMask[iy][ix] = NEGATIVE_PEDESTAL;
+                    continue; // NEGATIVE_PEDESTAL go to next pixel
-                    continue;
+                              // TODO! No pedestal update???
                 } else if (max > m_nSigma * rms) {
-                    eventMask[iy][ix] = PHOTON;
+                    event_type = eventType::PHOTON;
-
+                    if (value < max)
+                        continue; // Not max go to the next pixel
                 } else if (total > c3 * m_nSigma * rms) {
-                    eventMask[iy][ix] = PHOTON;
+                    event_type = eventType::PHOTON;
                 } else {
-                    if (late_update) {
+                    m_pedestal.push(iy, ix, frame(iy, ix));
-                        pedestal_updates.push_back({ix, iy, frame(iy, ix)});
+                    continue; // It was a pedestal value nothing to store
-                    } else {
-                        m_pedestal.push(iy, ix, frame(iy, ix));
-                    }
-                    continue;
                 }
-                if (eventMask[iy][ix] == PHOTON &&
-                    (frame(iy, ix) - m_pedestal.mean(iy, ix)) >= max) {
-                    eventMask[iy][ix] = PHOTON_MAX;
-                    DynamicCluster 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);
+                // Store cluster
-                         ir < (m_cluster_sizeY / 2) + 1; ir++) {
+                if (event_type == eventType::PHOTON && value >= max) {
-                        for (short ic = -(m_cluster_sizeX / 2);
+                    event_type = eventType::PHOTON_MAX;
-                             ic < (m_cluster_sizeX / 2) + 1; ic++) {
+
+                    short i = 0;
+                    std::vector<PEDESTAL_TYPE> cluster_data(m_cluster_sizeX *
+                                                            m_cluster_sizeY);
+
+                    for (short ir = -dy; ir < dy + 1; ir++) {
+                        for (short ic = -dx; ic < dx + 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)) -
                                     m_pedestal.mean(iy + ir, ix + ic);
-                                cluster.set<PEDESTAL_TYPE>(i, tmp);
+                                cluster_data[i] = tmp;
                                 i++;
                             }
                         }
                     }
-                    clusters.push_back(cluster);
+                    clusters.push_back(
+                        ix, iy,
+                        reinterpret_cast<std::byte *>(cluster_data.data()));
                 }
             }
         }
-        if (late_update) {
-            for (auto &update : pedestal_updates) {
-                m_pedestal.push(update.y, update.x, update.value);
-            }
-        }
         return clusters;
     }
 
@@ -176,7 +161,7 @@ class ClusterFinder {
         // 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;
+                eventMask[iy][ix] = eventType::PEDESTAL;
                 // initialize max and total
                 FRAME_TYPE max = std::numeric_limits<FRAME_TYPE>::min();
                 long double total = 0;
@@ -184,7 +169,7 @@ class ClusterFinder {
                     pedestal.push(iy, ix, frame(iy, ix));
                     continue;
                 }
-                eventMask[iy][ix] = NEIGHBOUR;
+                eventMask[iy][ix] = eventType::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++) {
@@ -220,18 +205,18 @@ class ClusterFinder {
                         tthr2 = tthr - tthr2;
                 }
                 if (total > tthr1 || max > tthr) {
-                    eventMask[iy][ix] = PHOTON;
+                    eventMask[iy][ix] = eventType::PHOTON;
                     nph(iy, ix) += 1;
                     rest(iy, ix) -= m_threshold;
                 } else {
                     pedestal.push(iy, ix, frame(iy, ix));
                     continue;
                 }
-                if (eventMask[iy][ix] == PHOTON &&
+                if (eventMask[iy][ix] == eventType::PHOTON &&
                     frame(iy, ix) - pedestal.mean(iy, ix) >= max) {
-                    eventMask[iy][ix] = PHOTON_MAX;
+                    eventMask[iy][ix] = eventType::PHOTON_MAX;
                     DynamicCluster cluster(m_cluster_sizeX, m_cluster_sizeY,
-                                    Dtype(typeid(FRAME_TYPE)));
+                                           Dtype(typeid(FRAME_TYPE)));
                     cluster.x = ix;
                     cluster.y = iy;
                     short i = 0;

include/aare/ClusterVector.hpp

@@ -0,0 +1,76 @@
+#pragma once
+#include <cstddef>
+#include <cstdint>
+
+template <typename T> class ClusterVector {
+    int32_t m_cluster_size_x;
+    int32_t m_cluster_size_y;
+    std::byte *m_data{};
+    size_t m_size{0};
+    size_t m_capacity{10};
+
+  public:
+    ClusterVector(int32_t cluster_size_x, int32_t cluster_size_y)
+        : m_cluster_size_x(cluster_size_x), m_cluster_size_y(cluster_size_y) {
+        size_t num_bytes = element_offset() * m_capacity;
+        m_data = new std::byte[num_bytes]{};
+        // fmt::print("Allocating {} bytes\n", num_bytes);
+    }
+
+    // data better hold data of the right size!
+    void push_back(int32_t x, int32_t y, const std::byte *data) {
+        if (m_size == m_capacity) {
+            m_capacity *= 2;
+            std::byte *new_data =
+                new std::byte[element_offset()*m_capacity]{};
+            std::copy(m_data,
+                      m_data + element_offset()*m_size,
+                      new_data);
+            delete[] m_data;
+            m_data = new_data;
+        }
+        std::byte *ptr = element_ptr(m_size);
+        *reinterpret_cast<int32_t *>(ptr) = x;
+        ptr += sizeof(int32_t);
+        *reinterpret_cast<int32_t *>(ptr) = y;
+        ptr += sizeof(int32_t);
+
+        std::copy(data, data + m_cluster_size_x * m_cluster_size_y * sizeof(T),
+                  ptr);
+        m_size++;
+    }
+
+    std::vector<T> sum(){
+        std::vector<T> sums(m_size);
+        for (size_t i = 0; i < m_size; i++) {
+            T sum = 0;
+            std::byte *ptr = element_ptr(i) + 2 * sizeof(int32_t);
+            for (size_t j = 0; j < m_cluster_size_x * m_cluster_size_y; j++) {
+                sum += *reinterpret_cast<T *>(ptr);
+                ptr += sizeof(T);
+            }
+            sums[i] = sum;
+        }
+        return sums;
+    }
+
+    size_t size() const { return m_size; }
+    size_t element_offset() const {
+        return sizeof(m_cluster_size_x) + sizeof(m_cluster_size_y) +
+               m_cluster_size_x * m_cluster_size_y * sizeof(T);
+    }
+    size_t element_offset(size_t i) const {
+        return element_offset() * i;
+    }
+
+    std::byte* element_ptr(size_t i) {
+        return m_data + element_offset(i);
+    }
+
+    int16_t cluster_size_x() const { return m_cluster_size_x; }
+    int16_t cluster_size_y() const { return m_cluster_size_y; }
+
+    std::byte *data() { return m_data; }
+
+    ~ClusterVector() { delete[] m_data; }
+};

include/aare/Pedestal.hpp

@@ -18,6 +18,8 @@ template <typename SUM_TYPE = double> class Pedestal {
 
     uint32_t m_samples;
     NDArray<uint32_t, 2> m_cur_samples;
+    
+    //TODO! in case of int needs to be changed to uint64_t
     NDArray<SUM_TYPE, 2> m_sum;
     NDArray<SUM_TYPE, 2> m_sum2;
 

include/aare/defs.hpp

@@ -47,7 +47,7 @@ class DynamicCluster {
     int cluster_sizeY;
     int16_t x;
     int16_t y;
-    Dtype dt;
+    Dtype dt; // 4 bytes
 
   private:
     std::byte *m_data;

python/examples/play.py

@@ -1,15 +1,53 @@
+import sys
+sys.path.append('/home/l_msdetect/erik/aare/build')
+
+#Our normal python imports
+from pathlib import Path
 import matplotlib.pyplot as plt
 import numpy as np
-plt.ion()
+import boost_histogram as bh
-from pathlib import Path
+import time
-from aare import ClusterFile
 
-base = Path('~/data/aare_test_data/clusters').expanduser()
+from aare import File, ClusterFinder, VarClusterFinder
 
-f = ClusterFile(base / 'beam_En700eV_-40deg_300V_10us_d0_f0_100.clust')
+base = Path('/mnt/sls_det_storage/matterhorn_data/aare_test_data/')
-# f = ClusterFile(base / 'single_frame_97_clustrers.clust')
+
+f = File(base/'Moench03new/cu_half_speed_master_4.json')
+cf = ClusterFinder((400,400), (3,3))
+for i in range(1000):
+    cf.push_pedestal_frame(f.read_frame())
+
+fig, ax = plt.subplots()
+im = ax.imshow(cf.pedestal())
+cf.pedestal()
+cf.noise()
+
+N = 200
+t0 = time.perf_counter()
+hist1 = bh.Histogram(bh.axis.Regular(40, -2, 4000))
+f.seek(0)
+
+t0 = time.perf_counter()
+data = f.read_n(N)
+t_elapsed = time.perf_counter()-t0
+
+print(f'Reading {N} frames took {t_elapsed:.3f}s {N/t_elapsed:.0f} FPS')
+
+clusters = []
+for frame in data:
+    clusters += [cf.find_clusters_without_threshold(frame)]
 
 
-for i in range(10):
+t_elapsed = time.perf_counter()-t0
-    fn, cl = f.read_frame()
+print(f'Clustering {N} frames took {t_elapsed:.2f}s  {N/t_elapsed:.0f} FPS')
-    print(fn, cl.size)
+
+
+t0 = time.perf_counter()
+total_clusters = 0
+for cl in clusters:
+        arr = np.array(cl, copy = False)
+        hist1.fill(arr['data'].sum(axis = 1).sum(axis = 1))
+        total_clusters += cl.size
+# t_elapsed = time.perf_counter()-t0
+print(f'Filling histogram with {total_clusters} clusters took: {t_elapsed:.3f}s')
+print(f'Cluster per frame {total_clusters/N:.3f}')

python/src/cluster.hpp

@@ -1,4 +1,5 @@
 #include "aare/ClusterFinder.hpp"
+#include "aare/ClusterVector.hpp"
 #include "aare/NDView.hpp"
 #include "aare/Pedestal.hpp"
 #include "np_helper.hpp"
@@ -9,30 +10,98 @@
 #include <pybind11/stl.h>
 
 namespace py = pybind11;
+using pd_type = double;
+
+template <typename T>
+void define_cluster_vector(py::module &m, const std::string &typestr) {
+    auto class_name = fmt::format("ClusterVector_{}", typestr);
+    py::class_<ClusterVector<T>>(m, class_name.c_str(), py::buffer_protocol())
+        .def(py::init<int, int>())
+        .def_property_readonly("size", &ClusterVector<T>::size)
+        .def("element_offset",
+             py::overload_cast<>(&ClusterVector<T>::element_offset, py::const_))
+        .def_property_readonly("fmt",
+             [typestr](ClusterVector<T> &v) {
+                 return fmt::format(
+                     "i:x:\ni:y:\n({},{}){}:data:", v.cluster_size_x(),
+                     v.cluster_size_y(), typestr);
+             })
+        .def("sum", &ClusterVector<T>::sum)
+        .def_buffer([typestr](ClusterVector<T> &v) -> py::buffer_info {
+            return py::buffer_info(
+                v.data(),           /* Pointer to buffer */
+                v.element_offset(), /* Size of one scalar */
+                fmt::format("i:x:\ni:y:\n({},{}){}:data:", v.cluster_size_x(),
+                            v.cluster_size_y(),
+                            typestr), /* Format descriptor */
+                1,                    /* Number of dimensions */
+                {v.size()},           /* Buffer dimensions */
+                {v.element_offset()}  /* Strides (in bytes) for each index */
+            );
+        });
+}
 
 void define_cluster_finder_bindings(py::module &m) {
-    py::class_<ClusterFinder<uint16_t, double>>(m, "ClusterFinder")
+    py::class_<ClusterFinder<uint16_t, pd_type>>(m, "ClusterFinder")
         .def(py::init<Shape<2>, Shape<2>>())
         .def("push_pedestal_frame",
-             [](ClusterFinder<uint16_t, double> &self,
+             [](ClusterFinder<uint16_t, pd_type> &self,
                 py::array_t<uint16_t> frame) {
                  auto view = make_view_2d(frame);
                  self.push_pedestal_frame(view);
              })
         .def("pedestal",
-             [](ClusterFinder<uint16_t, double> &self) {
+             [](ClusterFinder<uint16_t, pd_type> &self) {
-                 auto pd = new NDArray<double, 2>{};
+                 auto pd = new NDArray<pd_type, 2>{};
                  *pd = self.pedestal();
                  return return_image_data(pd);
              })
+        .def("noise",
+             [](ClusterFinder<uint16_t, pd_type> &self) {
+                 auto arr = new NDArray<pd_type, 2>{};
+                 *arr = self.noise();
+                 return return_image_data(arr);
+             })
         .def("find_clusters_without_threshold",
-             [](ClusterFinder<uint16_t, double> &self,
+             [](ClusterFinder<uint16_t, pd_type> &self,
                 py::array_t<uint16_t> frame) {
                  auto view = make_view_2d(frame);
-                 auto clusters = self.find_clusters_without_threshold(view);
+                 auto *vec = new ClusterVector<pd_type>(
-                 return clusters;
+                     self.find_clusters_without_threshold(view));
+                 return vec;
              });
 
+    m.def("hello", []() {
+        fmt::print("Hello from C++\n");
+        auto v = new ClusterVector<int>(3, 3);
+        int data[] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
+        v->push_back(5, 30, reinterpret_cast<std::byte *>(data));
+        v->push_back(5, 55, reinterpret_cast<std::byte *>(data));
+        v->push_back(5, 20, reinterpret_cast<std::byte *>(data));
+        v->push_back(5, 30, reinterpret_cast<std::byte *>(data));
+
+        return v;
+    });
+
+    define_cluster_vector<int>(m, "i");
+    define_cluster_vector<double>(m, "d");
+
+    // py::class_<ClusterVector<int>>(m, "ClusterVector", py::buffer_protocol())
+    //     .def(py::init<int, int>())
+    //     .def("size", &ClusterVector<int>::size)
+    //     .def("element_offset",
+    //     py::overload_cast<>(&ClusterVector<int>::element_offset, py::const_))
+    //     .def_buffer([](ClusterVector<int> &v) -> py::buffer_info {
+    //         return py::buffer_info(
+    //             v.data(),            /* Pointer to buffer */
+    //             v.element_offset(),  /* Size of one scalar */
+    //             fmt::format("h:x:\nh:y:\n({},{})i:data:", v.cluster_size_x(),
+    //             v.cluster_size_y()), /* Format descriptor */ 1, /* Number of
+    //             dimensions */ {v.size()},          /* Buffer dimensions */
+    //             {v.element_offset()} /* Strides (in bytes) for each index */
+    //         );
+    //     });
+
     py::class_<DynamicCluster>(m, "DynamicCluster", py::buffer_protocol())
         .def(py::init<int, int, Dtype>())
         .def("size", &DynamicCluster::size)