From 60534add92fefbc151cb0d4af571008d3db89426 Mon Sep 17 00:00:00 2001 From: froejdh_e Date: Wed, 11 Dec 2024 09:54:33 +0100 Subject: [PATCH] WIP --- include/aare/ClusterFile.hpp | 25 +++-- include/aare/ClusterVector.hpp | 46 ++++---- include/aare/File.hpp | 1 + python/src/cluster.hpp | 20 ++-- python/src/cluster_file.hpp | 33 ++++-- python/src/file.hpp | 3 +- src/ClusterFile.cpp | 192 +++++++++++++++++++-------------- src/File.cpp | 2 + 8 files changed, 198 insertions(+), 124 deletions(-) diff --git a/include/aare/ClusterFile.hpp b/include/aare/ClusterFile.hpp index f866dd6..edcb91e 100644 --- a/include/aare/ClusterFile.hpp +++ b/include/aare/ClusterFile.hpp @@ -1,6 +1,7 @@ #pragma once #include "aare/defs.hpp" +#include "aare/ClusterVector.hpp" #include #include @@ -38,30 +39,40 @@ struct ClusterAnalysis { double etay; }; - - +/* +Binary cluster file. Expects data to be layed out as: +int32_t frame_number +uint32_t number_of_clusters +int16_t x, int16_t y, int32_t data[9] x number_of_clusters +int32_t frame_number +uint32_t number_of_clusters +.... +*/ class ClusterFile { FILE *fp{}; uint32_t m_num_left{}; size_t m_chunk_size{}; + const std::string m_mode; public: - ClusterFile(const std::filesystem::path &fname, size_t chunk_size = 1000); + ClusterFile(const std::filesystem::path &fname, size_t chunk_size = 1000, + const std::string &mode = "r"); ~ClusterFile(); std::vector read_clusters(size_t n_clusters); std::vector read_frame(int32_t &out_fnum); + void write_frame(int32_t frame_number, const ClusterVector& clusters); std::vector read_cluster_with_cut(size_t n_clusters, double *noise_map, int nx, int ny); int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, - double *eta2x, double *eta2y, double *eta3x, double *eta3y); + double *eta2x, double *eta2y, double *eta3x, + double *eta3y); int analyze_cluster(Cluster cl, int32_t *t2, int32_t *t3, char *quad, - double *eta2x, double *eta2y, double *eta3x, - double *eta3y); + double *eta2x, double *eta2y, double *eta3x, + double *eta3y); size_t chunk_size() const { return m_chunk_size; } void close(); - }; } // namespace aare diff --git a/include/aare/ClusterVector.hpp b/include/aare/ClusterVector.hpp index 5445cbf..76e7e21 100644 --- a/include/aare/ClusterVector.hpp +++ b/include/aare/ClusterVector.hpp @@ -6,14 +6,25 @@ #include template class ClusterVector { - int32_t m_cluster_size_x; - int32_t m_cluster_size_y; + using value_type = T; + using coord_t = int16_t; + coord_t m_cluster_size_x; + coord_t m_cluster_size_y; std::byte *m_data{}; size_t m_size{0}; size_t m_capacity; + /* + 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: - ClusterVector(int32_t cluster_size_x, int32_t cluster_size_y, + ClusterVector(coord_t cluster_size_x, coord_t cluster_size_y, size_t capacity = 1024) : m_cluster_size_x(cluster_size_x), m_cluster_size_y(cluster_size_y), m_capacity(capacity) { @@ -25,15 +36,8 @@ template class ClusterVector { delete[] m_data; } - // ClusterVector(const ClusterVector & other){ - // m_cluster_size_x = other.m_cluster_size_x; - // m_cluster_size_y = other.m_cluster_size_y; - // m_data = new std::byte[other.m_capacity]; - // std::copy(other.m_data, other.m_data + other.m_capacity, m_data); - // m_size = other.m_size; - // m_capacity = other.m_capacity; - // } - + + //Move constructor ClusterVector(ClusterVector &&other) noexcept : m_cluster_size_x(other.m_cluster_size_x), m_cluster_size_y(other.m_cluster_size_y), m_data(other.m_data), @@ -66,15 +70,15 @@ template class ClusterVector { } // data better hold data of the right size! - void push_back(int32_t x, int32_t y, const std::byte *data) { + void push_back(coord_t x, coord_t y, const std::byte *data) { if (m_size == m_capacity) { allocate_buffer(m_capacity * 2); } std::byte *ptr = element_ptr(m_size); - *reinterpret_cast(ptr) = x; - ptr += sizeof(int32_t); - *reinterpret_cast(ptr) = y; - ptr += sizeof(int32_t); + *reinterpret_cast(ptr) = x; + ptr += sizeof(coord_t); + *reinterpret_cast(ptr) = y; + ptr += sizeof(coord_t); std::copy(data, data + m_cluster_size_x * m_cluster_size_y * sizeof(T), ptr); @@ -85,7 +89,7 @@ template class ClusterVector { std::vector sums(m_size); const size_t stride = element_offset(); const size_t n_pixels = m_cluster_size_x * m_cluster_size_y; - std::byte *ptr = m_data + 2 * sizeof(int32_t); // skip x and y + std::byte *ptr = m_data + 2 * sizeof(coord_t); // skip x and y for (size_t i = 0; i < m_size; i++) { sums[i] = @@ -109,6 +113,12 @@ template class ClusterVector { int16_t cluster_size_y() const { return m_cluster_size_y; } std::byte *data() { return m_data; } + const std::byte *data() const { return m_data; } + + const std::string_view fmt_base() const { + //TODO! how do we match on coord_t? + return m_fmt_base; + } private: void allocate_buffer(size_t new_capacity) { diff --git a/include/aare/File.hpp b/include/aare/File.hpp index b29171a..7aa30e1 100644 --- a/include/aare/File.hpp +++ b/include/aare/File.hpp @@ -44,6 +44,7 @@ class File { void read_into(std::byte *image_buf); void read_into(std::byte *image_buf, size_t n_frames); + size_t frame_number(); //!< get the frame number at the current position size_t frame_number(size_t frame_index); //!< get the frame number at the given frame index size_t bytes_per_frame() const; size_t pixels_per_frame() const; diff --git a/python/src/cluster.hpp b/python/src/cluster.hpp index 480aea1..0fef093 100644 --- a/python/src/cluster.hpp +++ b/python/src/cluster.hpp @@ -21,25 +21,25 @@ void define_cluster_vector(py::module &m, const std::string &typestr) { .def("element_offset", py::overload_cast<>(&ClusterVector::element_offset, py::const_)) .def_property_readonly("fmt", - [typestr](ClusterVector &v) { + [typestr](ClusterVector &self) { return fmt::format( - "i:x:\ni:y:\n({},{}){}:data:", v.cluster_size_x(), - v.cluster_size_y(), typestr); + self.fmt_base(), self.cluster_size_x(), + self.cluster_size_y(), typestr); }) .def("sum", [](ClusterVector &self) { auto *vec = new std::vector(self.sum()); return return_vector(vec); }) - .def_buffer([typestr](ClusterVector &v) -> py::buffer_info { + .def_buffer([typestr](ClusterVector &self) -> 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(), + self.data(), /* Pointer to buffer */ + self.element_offset(), /* Size of one scalar */ + fmt::format(self.fmt_base(), self.cluster_size_x(), + self.cluster_size_y(), typestr), /* Format descriptor */ 1, /* Number of dimensions */ - {v.size()}, /* Buffer dimensions */ - {v.element_offset()} /* Strides (in bytes) for each index */ + {self.size()}, /* Buffer dimensions */ + {self.element_offset()} /* Strides (in bytes) for each index */ ); }); } diff --git a/python/src/cluster_file.hpp b/python/src/cluster_file.hpp index 6f37c3d..543073f 100644 --- a/python/src/cluster_file.hpp +++ b/python/src/cluster_file.hpp @@ -1,7 +1,6 @@ #include "aare/ClusterFile.hpp" #include "aare/defs.hpp" - #include #include #include @@ -18,33 +17,47 @@ void define_cluster_file_io_bindings(py::module &m) { PYBIND11_NUMPY_DTYPE(Cluster, x, y, data); py::class_(m, "ClusterFile") - .def(py::init(), py::arg(), py::arg("chunk_size") = 1000) + .def(py::init(), + py::arg(), py::arg("chunk_size") = 1000, py::arg("mode") = "r") .def("read_clusters", [](ClusterFile &self, size_t n_clusters) { - auto* vec = new std::vector(self.read_clusters(n_clusters)); + auto *vec = + new std::vector(self.read_clusters(n_clusters)); return return_vector(vec); }) .def("read_frame", [](ClusterFile &self) { int32_t frame_number; - auto* vec = new std::vector(self.read_frame(frame_number)); + auto *vec = + new std::vector(self.read_frame(frame_number)); return py::make_tuple(frame_number, return_vector(vec)); }) + .def("write_frame", &ClusterFile::write_frame) .def("read_cluster_with_cut", - [](ClusterFile &self, size_t n_clusters, py::array_t noise_map, int nx, int ny) { + [](ClusterFile &self, size_t n_clusters, + py::array_t noise_map, int nx, int ny) { auto view = make_view_2d(noise_map); - auto* vec = new std::vector(self.read_cluster_with_cut(n_clusters, view.data(), nx, ny)); + auto *vec = + new std::vector(self.read_cluster_with_cut( + n_clusters, view.data(), nx, ny)); return return_vector(vec); }) .def("__enter__", [](ClusterFile &self) { return &self; }) - .def("__exit__", [](ClusterFile &self) { self.close();}) + .def("__exit__", + [](ClusterFile &self, + const std::optional &exc_type, + const std::optional &exc_value, + const std::optional &traceback) { + self.close(); + }) .def("__iter__", [](ClusterFile &self) { return &self; }) .def("__next__", [](ClusterFile &self) { - auto vec = new std::vector(self.read_clusters(self.chunk_size())); - if(vec->size() == 0) { + auto vec = + new std::vector(self.read_clusters(self.chunk_size())); + if (vec->size() == 0) { throw py::stop_iteration(); } return return_vector(vec); }); - } \ No newline at end of file diff --git a/python/src/file.hpp b/python/src/file.hpp index 3c44c43..30fa82f 100644 --- a/python/src/file.hpp +++ b/python/src/file.hpp @@ -51,7 +51,8 @@ void define_file_io_bindings(py::module &m) { .def(py::init()) - .def("frame_number", &File::frame_number) + .def("frame_number", py::overload_cast<>(&File::frame_number)) + .def("frame_number", py::overload_cast(&File::frame_number)) .def_property_readonly("bytes_per_frame", &File::bytes_per_frame) .def_property_readonly("pixels_per_frame", &File::pixels_per_frame) .def("seek", &File::seek) diff --git a/src/ClusterFile.cpp b/src/ClusterFile.cpp index 3daa9d6..182726b 100644 --- a/src/ClusterFile.cpp +++ b/src/ClusterFile.cpp @@ -2,25 +2,54 @@ namespace aare { -ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size): m_chunk_size(chunk_size) { - fp = fopen(fname.c_str(), "rb"); - if (!fp) { - throw std::runtime_error("Could not open file: " + fname.string()); +ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size, + const std::string &mode) + : m_chunk_size(chunk_size), m_mode(mode) { + + if (mode == "r") { + fp = fopen(fname.c_str(), "rb"); + if (!fp) { + throw std::runtime_error("Could not open file for reading: " + fname.string()); + } + } else if (mode == "w") { + fp = fopen(fname.c_str(), "wb"); + if (!fp) { + throw std::runtime_error("Could not open file for writing: " + fname.string()); + } + } else { + throw std::runtime_error("Unsupported mode: " + mode); + } + +} + +ClusterFile::~ClusterFile() { close(); } + +void ClusterFile::close() { + if (fp) { + fclose(fp); + fp = nullptr; } } -ClusterFile::~ClusterFile() { - close(); -} - -void ClusterFile::close(){ - if (fp){ - fclose(fp); - fp = nullptr; - } +void ClusterFile::write_frame(int32_t frame_number, const ClusterVector& clusters){ + if (m_mode != "w") { + 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"); + } + fwrite(&frame_number, sizeof(frame_number), 1, fp); + uint32_t n_clusters = clusters.size(); + fwrite(&n_clusters, sizeof(n_clusters), 1, fp); + fwrite(clusters.data(), clusters.element_offset(), clusters.size(), fp); + // write clusters + // fwrite(clusters.data(), sizeof(Cluster), clusters.size(), fp); } std::vector ClusterFile::read_clusters(size_t n_clusters) { + if (m_mode != "r") { + throw std::runtime_error("File not opened for reading"); + } std::vector clusters(n_clusters); int32_t iframe = 0; // frame number needs to be 4 bytes! @@ -38,7 +67,8 @@ std::vector ClusterFile::read_clusters(size_t n_clusters) { } else { nn = nph; } - nph_read += fread(reinterpret_cast(buf + nph_read), sizeof(Cluster), nn, fp); + nph_read += fread(reinterpret_cast(buf + nph_read), + sizeof(Cluster), nn, fp); m_num_left = nph - nn; // write back the number of photons left } @@ -52,8 +82,8 @@ std::vector ClusterFile::read_clusters(size_t n_clusters) { else nn = nph; - nph_read += - fread(reinterpret_cast(buf + nph_read), sizeof(Cluster), nn, fp); + nph_read += fread(reinterpret_cast(buf + nph_read), + sizeof(Cluster), nn, fp); m_num_left = nph - nn; } if (nph_read >= n_clusters) @@ -68,8 +98,12 @@ std::vector ClusterFile::read_clusters(size_t n_clusters) { } std::vector ClusterFile::read_frame(int32_t &out_fnum) { + if (m_mode != "r") { + throw std::runtime_error("File not opened for reading"); + } if (m_num_left) { - throw std::runtime_error("There are still photons left in the last frame"); + throw std::runtime_error( + "There are still photons left in the last frame"); } if (fread(&out_fnum, sizeof(out_fnum), 1, fp) != 1) { @@ -82,17 +116,19 @@ std::vector ClusterFile::read_frame(int32_t &out_fnum) { } std::vector clusters(n_clusters); - if (fread(clusters.data(), sizeof(Cluster), n_clusters, fp) != static_cast(n_clusters)) { + if (fread(clusters.data(), sizeof(Cluster), n_clusters, fp) != + static_cast(n_clusters)) { throw std::runtime_error("Could not read clusters"); } return clusters; - } std::vector ClusterFile::read_cluster_with_cut(size_t n_clusters, double *noise_map, int nx, int ny) { - + if (m_mode != "r") { + throw std::runtime_error("File not opened for reading"); + } std::vector clusters(n_clusters); // size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf, // uint32_t *n_left, double *noise_map, int @@ -124,7 +160,8 @@ std::vector ClusterFile::read_cluster_with_cut(size_t n_clusters, } for (size_t iph = 0; iph < nn; iph++) { // read photons 1 by 1 - size_t n_read = fread(reinterpret_cast(ptr), sizeof(Cluster), 1, fp); + size_t n_read = + fread(reinterpret_cast(ptr), sizeof(Cluster), 1, fp); if (n_read != 1) { clusters.resize(nph_read); return clusters; @@ -158,71 +195,71 @@ std::vector ClusterFile::read_cluster_with_cut(size_t n_clusters, break; } } - if (nph_read < n_clusters) { - // // keep on reading frames and photons until reaching n_clusters - while (fread(&iframe, sizeof(iframe), 1, fp)) { - // // printf("%d\n",nph_read); + if (nph_read < n_clusters) { + // // keep on reading frames and photons until reaching + // n_clusters + while (fread(&iframe, sizeof(iframe), 1, fp)) { + // // printf("%d\n",nph_read); - if (fread(&nph, sizeof(nph), 1, fp)) { - // // printf("** %d\n",nph); - m_num_left = nph; - for (size_t iph = 0; iph < nph; iph++) { - // // read photons 1 by 1 - size_t n_read = - fread(reinterpret_cast(ptr), sizeof(Cluster), 1, fp); - if (n_read != 1) { - clusters.resize(nph_read); - return clusters; - // return nph_read; - } - good = 1; - if (noise_map) { - if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 && - ptr->y < ny) { - tot1 = ptr->data[4]; - analyze_cluster(*ptr, &t2max, &tot3, NULL, - NULL, - NULL, NULL, NULL); - // noise = noise_map[ptr->y * nx + ptr->x]; - noise = noise_map[ptr->y + ny * ptr->x]; - if (tot1 > noise || t2max > 2 * noise || - tot3 > 3 * noise) { - ; - } else - good = 0; - } else { - printf("Bad pixel number %d %d\n", ptr->x, - ptr->y); good = 0; - } - } - if (good) { - ptr++; - nph_read++; - } - (m_num_left)--; - if (nph_read >= n_clusters) - break; + if (fread(&nph, sizeof(nph), 1, fp)) { + // // printf("** %d\n",nph); + m_num_left = nph; + for (size_t iph = 0; iph < nph; iph++) { + // // read photons 1 by 1 + size_t n_read = fread(reinterpret_cast(ptr), + sizeof(Cluster), 1, fp); + if (n_read != 1) { + clusters.resize(nph_read); + return clusters; + // return nph_read; } + good = 1; + if (noise_map) { + if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 && + ptr->y < ny) { + tot1 = ptr->data[4]; + analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL, + NULL, NULL, NULL); + // noise = noise_map[ptr->y * nx + ptr->x]; + noise = noise_map[ptr->y + ny * ptr->x]; + if (tot1 > noise || t2max > 2 * noise || + tot3 > 3 * noise) { + ; + } else + good = 0; + } else { + printf("Bad pixel number %d %d\n", ptr->x, ptr->y); + good = 0; + } + } + if (good) { + ptr++; + nph_read++; + } + (m_num_left)--; + if (nph_read >= n_clusters) + break; } - if (nph_read >= n_clusters) - break; } + if (nph_read >= n_clusters) + break; } - // printf("%d\n",nph_read); - clusters.resize(nph_read); - return clusters; - + } + // printf("%d\n",nph_read); + clusters.resize(nph_read); + return clusters; } -int ClusterFile::analyze_cluster(Cluster cl, int32_t *t2, int32_t *t3, char *quad, - double *eta2x, double *eta2y, double *eta3x, - double *eta3y) { +int ClusterFile::analyze_cluster(Cluster cl, int32_t *t2, int32_t *t3, + char *quad, double *eta2x, double *eta2y, + double *eta3x, double *eta3y) { return analyze_data(cl.data, t2, t3, quad, eta2x, eta2y, eta3x, eta3y); } -int ClusterFile::analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, - double *eta2x, double *eta2y, double *eta3x, double *eta3y) { +int ClusterFile::analyze_data(int32_t *data, int32_t *t2, int32_t *t3, + char *quad, double *eta2x, double *eta2y, + double *eta3x, double *eta3y) { int ok = 1; @@ -263,7 +300,8 @@ int ClusterFile::analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *qua c = i; } } - //printf("*** %d %d %d %d -- %d\n",tot2[0],tot2[1],tot2[2],tot2[3],t2max); + // printf("*** %d %d %d %d -- + // %d\n",tot2[0],tot2[1],tot2[2],tot2[3],t2max); if (quad) *quad = c; if (t2) @@ -318,6 +356,4 @@ int ClusterFile::analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *qua return ok; } - - } // namespace aare \ No newline at end of file diff --git a/src/File.cpp b/src/File.cpp index d45e903..37e4c57 100644 --- a/src/File.cpp +++ b/src/File.cpp @@ -58,6 +58,8 @@ void File::read_into(std::byte *image_buf) { file_impl->read_into(image_buf); } void File::read_into(std::byte *image_buf, size_t n_frames) { file_impl->read_into(image_buf, n_frames); } + +size_t File::frame_number() { return file_impl->frame_number(tell()); } size_t File::frame_number(size_t frame_index) { return file_impl->frame_number(frame_index); }