From 6a150e8d98deaf5e29bf2cbd7d3ef485f5761616 Mon Sep 17 00:00:00 2001 From: froejdh_e Date: Tue, 10 Dec 2024 17:21:05 +0100 Subject: [PATCH] WIP --- include/aare/ClusterFile.hpp | 2 +- include/aare/ClusterFinder.hpp | 145 +++++++++++++++------------------ include/aare/ClusterVector.hpp | 76 +++++++++++++++++ include/aare/Pedestal.hpp | 2 + include/aare/defs.hpp | 2 +- python/examples/play.py | 56 +++++++++++-- python/src/cluster.hpp | 83 +++++++++++++++++-- 7 files changed, 268 insertions(+), 98 deletions(-) create mode 100644 include/aare/ClusterVector.hpp diff --git a/include/aare/ClusterFile.hpp b/include/aare/ClusterFile.hpp index 2baf0f4..f866dd6 100644 --- a/include/aare/ClusterFile.hpp +++ b/include/aare/ClusterFile.hpp @@ -1,4 +1,4 @@ - +#pragma once #include "aare/defs.hpp" #include diff --git a/include/aare/ClusterFinder.hpp b/include/aare/ClusterFinder.hpp index addb6db..33a00ea 100644 --- a/include/aare/ClusterFinder.hpp +++ b/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 m_pedestal; public: - ClusterFinder(Shape<2> image_size, Shape<2>cluster_size, 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_threshold(threshold), m_nSigma(nSigma), + ClusterFinder(Shape<2> image_size, Shape<2> cluster_size, + 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_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]) { - - // c2 = sqrt((cluster_sizeY + 1) / 2 * (cluster_sizeX + 1) / 2); - // c3 = sqrt(cluster_sizeX * cluster_sizeY); - }; + fmt::print("TypeIndex: {}\n", sizeof(Dtype)); + }; void push_pedestal_frame(NDView frame) { m_pedestal.push(frame); } - NDArray pedestal() { - return m_pedestal.mean(); - } + NDArray pedestal() { return m_pedestal.mean(); } - std::vector - find_clusters_without_threshold(NDView frame, - // Pedestal &pedestal, - bool late_update = false) { - struct pedestal_update { - int x; - int y; - FRAME_TYPE value; - }; - std::vector pedestal_updates; + NDArray noise() { return m_pedestal.std(); } - std::vector clusters; - std::vector> eventMask; - for (int i = 0; i < frame.shape(0); i++) { - eventMask.push_back(std::vector(frame.shape(1))); - } - long double val; - long double max; + ClusterVector + find_clusters_without_threshold(NDView frame) { + // std::vector clusters; + // std::vector clusters; //Hard coded 3x3 cluster + // clusters.reserve(2000); + ClusterVector clusters(m_cluster_sizeX, m_cluster_sizeY); + 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 - max = std::numeric_limits::min(); - long double total = 0; - eventMask[iy][ix] = PEDESTAL; + PEDESTAL_TYPE max = std::numeric_limits::min(); + PEDESTAL_TYPE total = 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++) { + 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)) { - val = frame(iy + ir, ix + ic) - - m_pedestal.mean(iy + ir, ix + ic); + PEDESTAL_TYPE val = + frame(iy + ir, ix + ic) - + m_pedestal.mean(iy + ir, ix + ic); + total += val; - if (val > max) { - max = val; - } + max = std::max(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) { - eventMask[iy][ix] = NEGATIVE_PEDESTAL; - continue; + if (value < -m_nSigma * rms) { + continue; // NEGATIVE_PEDESTAL go to next pixel + // 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) { - pedestal_updates.push_back({ix, iy, frame(iy, ix)}); - } else { - m_pedestal.push(iy, ix, frame(iy, ix)); - } - continue; + m_pedestal.push(iy, ix, frame(iy, ix)); + continue; // It was a pedestal value nothing to store } - 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); - ir < (m_cluster_sizeY / 2) + 1; ir++) { - for (short ic = -(m_cluster_sizeX / 2); - ic < (m_cluster_sizeX / 2) + 1; ic++) { + // Store cluster + if (event_type == eventType::PHOTON && value >= max) { + event_type = eventType::PHOTON_MAX; + + short i = 0; + std::vector 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( frame(iy + ir, ix + ic)) - m_pedestal.mean(iy + ir, ix + ic); - cluster.set(i, tmp); + cluster_data[i] = tmp; i++; } } } - clusters.push_back(cluster); + clusters.push_back( + ix, iy, + reinterpret_cast(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::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; diff --git a/include/aare/ClusterVector.hpp b/include/aare/ClusterVector.hpp new file mode 100644 index 0000000..6016e87 --- /dev/null +++ b/include/aare/ClusterVector.hpp @@ -0,0 +1,76 @@ +#pragma once +#include +#include + +template 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(ptr) = x; + ptr += sizeof(int32_t); + *reinterpret_cast(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 sum(){ + std::vector 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(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; } +}; \ No newline at end of file diff --git a/include/aare/Pedestal.hpp b/include/aare/Pedestal.hpp index b5f245b..216c204 100644 --- a/include/aare/Pedestal.hpp +++ b/include/aare/Pedestal.hpp @@ -18,6 +18,8 @@ template class Pedestal { uint32_t m_samples; NDArray m_cur_samples; + + //TODO! in case of int needs to be changed to uint64_t NDArray m_sum; NDArray m_sum2; diff --git a/include/aare/defs.hpp b/include/aare/defs.hpp index 35b9624..7466410 100644 --- a/include/aare/defs.hpp +++ b/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; diff --git a/python/examples/play.py b/python/examples/play.py index 633b7e2..c5abcbf 100644 --- a/python/examples/play.py +++ b/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() -from pathlib import Path -from aare import ClusterFile +import boost_histogram as bh +import time -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') -# f = ClusterFile(base / 'single_frame_97_clustrers.clust') +base = Path('/mnt/sls_det_storage/matterhorn_data/aare_test_data/') + +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): - fn, cl = f.read_frame() - print(fn, cl.size) +t_elapsed = time.perf_counter()-t0 +print(f'Clustering {N} frames took {t_elapsed:.2f}s {N/t_elapsed:.0f} FPS') + + +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}') \ No newline at end of file diff --git a/python/src/cluster.hpp b/python/src/cluster.hpp index 6932281..b9296fd 100644 --- a/python/src/cluster.hpp +++ b/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 namespace py = pybind11; +using pd_type = double; + +template +void define_cluster_vector(py::module &m, const std::string &typestr) { + auto class_name = fmt::format("ClusterVector_{}", typestr); + py::class_>(m, class_name.c_str(), py::buffer_protocol()) + .def(py::init()) + .def_property_readonly("size", &ClusterVector::size) + .def("element_offset", + py::overload_cast<>(&ClusterVector::element_offset, py::const_)) + .def_property_readonly("fmt", + [typestr](ClusterVector &v) { + return fmt::format( + "i:x:\ni:y:\n({},{}){}:data:", v.cluster_size_x(), + v.cluster_size_y(), typestr); + }) + .def("sum", &ClusterVector::sum) + .def_buffer([typestr](ClusterVector &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_>(m, "ClusterFinder") + py::class_>(m, "ClusterFinder") .def(py::init, Shape<2>>()) .def("push_pedestal_frame", - [](ClusterFinder &self, + [](ClusterFinder &self, py::array_t frame) { auto view = make_view_2d(frame); self.push_pedestal_frame(view); }) .def("pedestal", - [](ClusterFinder &self) { - auto pd = new NDArray{}; + [](ClusterFinder &self) { + auto pd = new NDArray{}; *pd = self.pedestal(); return return_image_data(pd); }) + .def("noise", + [](ClusterFinder &self) { + auto arr = new NDArray{}; + *arr = self.noise(); + return return_image_data(arr); + }) .def("find_clusters_without_threshold", - [](ClusterFinder &self, + [](ClusterFinder &self, py::array_t frame) { auto view = make_view_2d(frame); - auto clusters = self.find_clusters_without_threshold(view); - return clusters; + auto *vec = new ClusterVector( + self.find_clusters_without_threshold(view)); + return vec; }); + m.def("hello", []() { + fmt::print("Hello from C++\n"); + auto v = new ClusterVector(3, 3); + int data[] = {1, 2, 3, 4, 5, 6, 7, 8, 9}; + v->push_back(5, 30, reinterpret_cast(data)); + v->push_back(5, 55, reinterpret_cast(data)); + v->push_back(5, 20, reinterpret_cast(data)); + v->push_back(5, 30, reinterpret_cast(data)); + + return v; + }); + + define_cluster_vector(m, "i"); + define_cluster_vector(m, "d"); + + // py::class_>(m, "ClusterVector", py::buffer_protocol()) + // .def(py::init()) + // .def("size", &ClusterVector::size) + // .def("element_offset", + // py::overload_cast<>(&ClusterVector::element_offset, py::const_)) + // .def_buffer([](ClusterVector &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_(m, "DynamicCluster", py::buffer_protocol()) .def(py::init()) .def("size", &DynamicCluster::size)