From 6e7e81b36ba7ab9276405df004c792d275111fe8 Mon Sep 17 00:00:00 2001 From: AliceMazzoleni99 Date: Fri, 21 Mar 2025 16:32:54 +0100 Subject: [PATCH] complete mess but need to install RedHat 9 --- include/aare/ClusterFile.hpp | 48 +++++---- include/aare/ClusterVector.hpp | 48 +++++---- src/ClusterFile.cpp | 184 ++++++++++++++++++++------------- 3 files changed, 162 insertions(+), 118 deletions(-) diff --git a/include/aare/ClusterFile.hpp b/include/aare/ClusterFile.hpp index 5bea342..d35f362 100644 --- a/include/aare/ClusterFile.hpp +++ b/include/aare/ClusterFile.hpp @@ -8,16 +8,12 @@ namespace aare { -//TODO! Template this? -struct Cluster3x3 { - int16_t x; - int16_t y; - int32_t data[9]; -}; -struct Cluster2x2 { - int16_t x; - int16_t y; - int32_t data[4]; +template +struct Cluster { + CoordType x; + CoordType y; + T data[ClusterSizeX * ClusterSizeY]; }; typedef enum { @@ -93,8 +89,7 @@ class ClusterFile { */ ClusterFile(const std::filesystem::path &fname, size_t chunk_size = 1000, const std::string &mode = "r"); - - + ~ClusterFile(); /** @@ -109,26 +104,26 @@ class ClusterFile { /** * @brief Read a single frame from the file and return the clusters. The * cluster vector will have the frame number set. - * @throws std::runtime_error if the file is not opened for reading or the file pointer not - * at the beginning of a frame + * @throws std::runtime_error if the file is not opened for reading or the + * file pointer not at the beginning of a frame */ ClusterVector read_frame(); - void write_frame(const ClusterVector &clusters); - + // Need to be migrated to support NDArray and return a ClusterVector // std::vector - // read_cluster_with_cut(size_t n_clusters, double *noise_map, int nx, int ny); + // read_cluster_with_cut(size_t n_clusters, double *noise_map, int nx, int + // ny); /** * @brief Return the chunk size */ size_t chunk_size() const { return m_chunk_size; } - - + /** - * @brief Close the file. If not closed the file will be closed in the destructor + * @brief Close the file. If not closed the file will be closed in the + * destructor */ void close(); }; @@ -138,8 +133,17 @@ int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, int analyze_cluster(Cluster3x3 &cl, int32_t *t2, int32_t *t3, char *quad, double *eta2x, double *eta2y, double *eta3x, double *eta3y); -NDArray calculate_eta2(ClusterVector &clusters); -Eta2 calculate_eta2(Cluster3x3 &cl); +template +NDArray calculate_eta2(ClusterVector &clusters); + +template Eta2 calculate_eta2(Cluster &cl); + Eta2 calculate_eta2(Cluster2x2 &cl); +template Eta2 calculate_eta2(ClusterType &cl); + +template +Eta2 calculate_eta2(Cluster &cl); + } // namespace aare diff --git a/include/aare/ClusterVector.hpp b/include/aare/ClusterVector.hpp index 1c15a22..c5e66b7 100644 --- a/include/aare/ClusterVector.hpp +++ b/include/aare/ClusterVector.hpp @@ -10,6 +10,8 @@ namespace aare { +template 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 @@ -21,10 +23,12 @@ namespace aare { * @tparam CoordType data type of the x and y coordinates of the cluster * (normally int16_t) */ -template class ClusterVector { +template +class ClusterVector> { using value_type = T; - size_t m_cluster_size_x; - size_t m_cluster_size_y; + // size_t m_cluster_size_x; + // size_t m_cluster_size_y; std::byte *m_data{}; size_t m_size{0}; size_t m_capacity; @@ -40,6 +44,8 @@ template class ClusterVector { constexpr static char m_fmt_base[] = "=h:x:\nh:y:\n({},{}){}:data:"; public: + using ClusterType = Cluster; + /** * @brief Construct a new ClusterVector object * @param cluster_size_x size of the cluster in x direction @@ -48,10 +54,8 @@ template class ClusterVector { * @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 cluster_size_x = 3, size_t cluster_size_y = 3, - size_t capacity = 1024, uint64_t frame_number = 0) - : m_cluster_size_x(cluster_size_x), m_cluster_size_y(cluster_size_y), - m_capacity(capacity), m_frame_number(frame_number) { + ClusterVector(size_t capacity = 1024, uint64_t frame_number = 0) + : m_capacity(capacity), m_frame_number(frame_number) { allocate_buffer(capacity); } @@ -59,10 +63,8 @@ template class ClusterVector { // 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), - m_size(other.m_size), m_capacity(other.m_capacity), - m_frame_number(other.m_frame_number) { + : 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; @@ -72,8 +74,6 @@ template class ClusterVector { ClusterVector &operator=(ClusterVector &&other) noexcept { if (this != &other) { delete[] m_data; - m_cluster_size_x = other.m_cluster_size_x; - m_cluster_size_y = other.m_cluster_size_y; m_data = other.m_data; m_size = other.m_size; m_capacity = other.m_capacity; @@ -116,8 +116,7 @@ template class ClusterVector { *reinterpret_cast(ptr) = y; ptr += sizeof(CoordType); - std::copy(data, data + m_cluster_size_x * m_cluster_size_y * sizeof(T), - ptr); + std::copy(data, data + ClusterSizeX * ClusterSizeY * sizeof(T), ptr); m_size++; } ClusterVector &operator+=(const ClusterVector &other) { @@ -137,7 +136,7 @@ template class ClusterVector { std::vector sum() { std::vector sums(m_size); const size_t stride = item_size(); - const size_t n_pixels = m_cluster_size_x * m_cluster_size_y; + 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++) { @@ -159,7 +158,7 @@ template class ClusterVector { std::vector sums(m_size); const size_t stride = item_size(); - if (m_cluster_size_x != 3 || m_cluster_size_y != 3) { + if (ClusterSizeX != 3 || ClusterSizeY != 3) { throw std::runtime_error( "Only 3x3 clusters are supported for the 2x2 sum."); } @@ -196,8 +195,7 @@ template class ClusterVector { * @brief Return the size in bytes of a single cluster */ size_t item_size() const { - return 2 * sizeof(CoordType) + - m_cluster_size_x * m_cluster_size_y * sizeof(T); + return 2 * sizeof(CoordType) + ClusterSizeX * ClusterSizeY * sizeof(T); } /** @@ -217,8 +215,8 @@ template class ClusterVector { return m_data + element_offset(i); } - size_t cluster_size_x() const { return m_cluster_size_x; } - size_t cluster_size_y() const { return m_cluster_size_y; } + // size_t cluster_size_x() const { return m_cluster_size_x; } + // size_t cluster_size_y() const { return m_cluster_size_y; } std::byte *data() { return m_data; } std::byte const *data() const { return m_data; } @@ -227,12 +225,12 @@ template class ClusterVector { * @brief Return a reference to the i-th cluster casted to type V * @tparam V type of the cluster */ - template V &at(size_t i) { - return *reinterpret_cast(element_ptr(i)); + ClusterType &at(size_t i) { + return *reinterpret_cast(element_ptr(i)); } - template const V &at(size_t i) const { - return *reinterpret_cast(element_ptr(i)); + const ClusterType &at(size_t i) const { + return *reinterpret_cast(element_ptr(i)); } const std::string_view fmt_base() const { diff --git a/src/ClusterFile.cpp b/src/ClusterFile.cpp index be3f607..c6ae470 100644 --- a/src/ClusterFile.cpp +++ b/src/ClusterFile.cpp @@ -59,8 +59,8 @@ ClusterVector ClusterFile::read_clusters(size_t n_clusters) { if (m_mode != "r") { throw std::runtime_error("File not opened for reading"); } - - ClusterVector clusters(3,3, n_clusters); + + ClusterVector clusters(3, 3, n_clusters); int32_t iframe = 0; // frame number needs to be 4 bytes! size_t nph_read = 0; @@ -78,7 +78,7 @@ ClusterVector ClusterFile::read_clusters(size_t n_clusters) { } else { nn = nph; } - nph_read += fread((buf + nph_read*clusters.item_size()), + nph_read += fread((buf + nph_read * clusters.item_size()), clusters.item_size(), nn, fp); m_num_left = nph - nn; // write back the number of photons left } @@ -93,7 +93,7 @@ ClusterVector ClusterFile::read_clusters(size_t n_clusters) { else nn = nph; - nph_read += fread((buf + nph_read*clusters.item_size()), + nph_read += fread((buf + nph_read * clusters.item_size()), clusters.item_size(), nn, fp); m_num_left = nph - nn; } @@ -112,8 +112,8 @@ ClusterVector ClusterFile::read_clusters(size_t n_clusters, ROI roi) { if (m_mode != "r") { throw std::runtime_error("File not opened for reading"); } - - ClusterVector clusters(3,3); + + ClusterVector clusters(3, 3); clusters.reserve(n_clusters); int32_t iframe = 0; // frame number needs to be 4 bytes! @@ -124,7 +124,7 @@ ClusterVector ClusterFile::read_clusters(size_t n_clusters, ROI roi) { // auto buf = reinterpret_cast(clusters.data()); // auto buf = clusters.data(); - Cluster3x3 tmp; //this would break if the cluster size changes + Cluster3x3 tmp; // this would break if the cluster size changes // if there are photons left from previous frame read them first if (nph) { @@ -135,13 +135,15 @@ ClusterVector ClusterFile::read_clusters(size_t n_clusters, ROI roi) { } else { nn = nph; } - //Read one cluster, in the ROI push back - // nph_read += fread((buf + nph_read*clusters.item_size()), - // clusters.item_size(), nn, fp); - for(size_t i = 0; i < nn; i++){ + // Read one cluster, in the ROI push back + // nph_read += fread((buf + nph_read*clusters.item_size()), + // clusters.item_size(), nn, fp); + for (size_t i = 0; i < nn; i++) { fread(&tmp, sizeof(tmp), 1, fp); - if(tmp.x >= roi.xmin && tmp.x <= roi.xmax && tmp.y >= roi.ymin && tmp.y <= roi.ymax){ - clusters.push_back(tmp.x, tmp.y, reinterpret_cast(tmp.data)); + if (tmp.x >= roi.xmin && tmp.x <= roi.xmax && tmp.y >= roi.ymin && + tmp.y <= roi.ymax) { + clusters.push_back(tmp.x, tmp.y, + reinterpret_cast(tmp.data)); nph_read++; } } @@ -161,10 +163,13 @@ ClusterVector ClusterFile::read_clusters(size_t n_clusters, ROI roi) { // nph_read += fread((buf + nph_read*clusters.item_size()), // clusters.item_size(), nn, fp); - for(size_t i = 0; i < nn; i++){ + for (size_t i = 0; i < nn; i++) { fread(&tmp, sizeof(tmp), 1, fp); - if(tmp.x >= roi.xmin && tmp.x <= roi.xmax && tmp.y >= roi.ymin && tmp.y <= roi.ymax){ - clusters.push_back(tmp.x, tmp.y, reinterpret_cast(tmp.data)); + if (tmp.x >= roi.xmin && tmp.x <= roi.xmax && + tmp.y >= roi.ymin && tmp.y <= roi.ymax) { + clusters.push_back( + tmp.x, tmp.y, + reinterpret_cast(tmp.data)); nph_read++; } } @@ -210,7 +215,6 @@ ClusterVector ClusterFile::read_frame() { return clusters; } - // std::vector ClusterFile::read_cluster_with_cut(size_t n_clusters, // double *noise_map, // int nx, int ny) { @@ -218,7 +222,8 @@ ClusterVector ClusterFile::read_frame() { // 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, +// // size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster +// *buf, // // uint32_t *n_left, double *noise_map, int // // nx, int ny) { // int iframe = 0; @@ -249,7 +254,8 @@ ClusterVector ClusterFile::read_frame() { // for (size_t iph = 0; iph < nn; iph++) { // // read photons 1 by 1 // size_t n_read = -// fread(reinterpret_cast(ptr), sizeof(Cluster3x3), 1, fp); +// fread(reinterpret_cast(ptr), sizeof(Cluster3x3), 1, +// fp); // if (n_read != 1) { // clusters.resize(nph_read); // return clusters; @@ -257,12 +263,15 @@ ClusterVector ClusterFile::read_frame() { // // TODO! error handling on read // good = 1; // if (noise_map) { -// if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 && ptr->y < ny) { +// 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, +// analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL, NULL, +// NULL, // NULL); // noise = noise_map[ptr->y * nx + ptr->x]; -// if (tot1 > noise || t2max > 2 * noise || tot3 > 3 * noise) { +// if (tot1 > noise || t2max > 2 * noise || tot3 > 3 * +// noise) { // ; // } else { // good = 0; @@ -316,8 +325,8 @@ ClusterVector ClusterFile::read_frame() { // } else // good = 0; // } else { -// printf("Bad pixel number %d %d\n", ptr->x, ptr->y); -// good = 0; +// printf("Bad pixel number %d %d\n", ptr->x, +// ptr->y); good = 0; // } // } // if (good) { @@ -338,37 +347,81 @@ ClusterVector ClusterFile::read_frame() { // return clusters; // } -NDArray calculate_eta2(ClusterVector &clusters) { - //TOTO! make work with 2x2 clusters +template +NDArray calculate_eta2(ClusterVector &clusters) { + // TOTO! make work with 2x2 clusters NDArray eta2({static_cast(clusters.size()), 2}); - - if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) { - for (size_t i = 0; i < clusters.size(); i++) { - auto e = calculate_eta2(clusters.at(i)); - eta2(i, 0) = e.x; - eta2(i, 1) = e.y; - } - }else if(clusters.cluster_size_x() == 2 || clusters.cluster_size_y() == 2){ - for (size_t i = 0; i < clusters.size(); i++) { - auto e = calculate_eta2(clusters.at(i)); - eta2(i, 0) = e.x; - eta2(i, 1) = e.y; - } - }else{ - throw std::runtime_error("Only 3x3 and 2x2 clusters are supported"); + + for (size_t i = 0; i < clusters.size(); i++) { + auto e = calculate_eta2(clusters.at(i)); + eta2(i, 0) = e.x; + eta2(i, 1) = e.y; } - + return eta2; } -/** - * @brief Calculate the eta2 values for a 3x3 cluster and return them in a Eta2 struct - * containing etay, etax and the corner of the cluster. -*/ -Eta2 calculate_eta2(Cluster3x3 &cl) { +/** + * @brief Calculate the eta2 values for a generic sized cluster and return them + * in a Eta2 struct containing etay, etax and the index of the respective 2x2 + * subcluster. + */ +template +Eta2 calculate_eta2(Cluster &cl) { Eta2 eta{}; - std::array tot2; + // TODO loads of overhead for a 2x2 clsuter maybe keep 2x2 calculation + size_t num_2x2_subclusters = (ClusterSizeX - 1) * (ClusterSizeY - 1); + std::array sum_2x2_subcluster; + for (size_t i = 0; i < ClusterSizeY - 1; ++i) { + for (size_t j = 0; j < ClusterSizeX - 1; ++j) + sum_2x2_subcluster[i * (ClusterSizeX - 1) + j] = + cl.data[i * ClusterSizeX + j] + + cl.data[i * ClusterSizeX + j + 1] + + cl.data[(i + 1) * ClusterSizeX + j] + + cl.data[(i + 1) * ClusterSizeX + j + 1]; + } + + auto c = std::max_element(sum_2x2_subclusters.begin(), + sum_2x2_subcluster.end()) - + sum_2x2_subcluster.begin(); + + eta.sum = sum_2x2_subcluster[c]; + + eta.x = static_cast(cl.data[(c + 1) * ClusterSizeX + 1]) / + (cl.data[0] + cl.data[1]); + + size_t index_top_left_2x2_subcluster = + (int(c / (ClusterSizeX - 1)) + 1) * ClusterSizeX + + c % (ClusterSizeX - 1) * 2 + 1; + if ((cl.data[index_top_left_2x2_subcluster] + + cl.data[index_top_left_2x2_subcluster - 1]) != 0) + eta.x = + static_cast(cl.data[index_top_left_2x2_subcluster] / + (cl.data[index_top_left_2x2_subcluster] + + cl.data[index_top_left_2x2_subcluster - 1])); + + if ((cl.data[index_top_left_2x2_subcluster] + + cl.data[index_top_left_2x2_subcluster - ClusterSizeX]) != 0) + eta.y = static_cast( + cl.data[index_top_left_2x2_subcluster] / + (cl.data[index_top_left_2x2_subcluster] + + cl.data[index_top_left_2x2_subcluster - ClusterSizeX])); + + eta.c = c; // TODO only supported for 2x2 and 3x3 clusters -> at least no + // underyling enum class + return eta; +} + +/** + * @brief Calculate the eta2 values for a 3x3 cluster and return them in a Eta2 + * struct containing etay, etax and the corner of the cluster. + */ +template Eta2 calculate_eta2(Cluster &cl) { + Eta2 eta{}; + + std::array tot2; tot2[0] = cl.data[0] + cl.data[1] + cl.data[3] + cl.data[4]; tot2[1] = cl.data[1] + cl.data[2] + cl.data[4] + cl.data[5]; tot2[2] = cl.data[3] + cl.data[4] + cl.data[6] + cl.data[7]; @@ -379,58 +432,47 @@ Eta2 calculate_eta2(Cluster3x3 &cl) { switch (c) { case cBottomLeft: if ((cl.data[3] + cl.data[4]) != 0) - eta.x = - static_cast(cl.data[4]) / (cl.data[3] + cl.data[4]); + eta.x = static_cast(cl.data[4]) / (cl.data[3] + cl.data[4]); if ((cl.data[1] + cl.data[4]) != 0) - eta.y = - static_cast(cl.data[4]) / (cl.data[1] + cl.data[4]); + eta.y = static_cast(cl.data[4]) / (cl.data[1] + cl.data[4]); eta.c = cBottomLeft; break; case cBottomRight: if ((cl.data[2] + cl.data[5]) != 0) - eta.x = - static_cast(cl.data[5]) / (cl.data[4] + cl.data[5]); + eta.x = static_cast(cl.data[5]) / (cl.data[4] + cl.data[5]); if ((cl.data[1] + cl.data[4]) != 0) - eta.y = - static_cast(cl.data[4]) / (cl.data[1] + cl.data[4]); + eta.y = static_cast(cl.data[4]) / (cl.data[1] + cl.data[4]); eta.c = cBottomRight; break; case cTopLeft: if ((cl.data[7] + cl.data[4]) != 0) - eta.x = - static_cast(cl.data[4]) / (cl.data[3] + cl.data[4]); + eta.x = static_cast(cl.data[4]) / (cl.data[3] + cl.data[4]); if ((cl.data[7] + cl.data[4]) != 0) - eta.y = - static_cast(cl.data[7]) / (cl.data[7] + cl.data[4]); + eta.y = static_cast(cl.data[7]) / (cl.data[7] + cl.data[4]); eta.c = cTopLeft; break; case cTopRight: if ((cl.data[5] + cl.data[4]) != 0) - eta.x = - static_cast(cl.data[5]) / (cl.data[5] + cl.data[4]); + eta.x = static_cast(cl.data[5]) / (cl.data[5] + cl.data[4]); if ((cl.data[7] + cl.data[4]) != 0) - eta.y = - static_cast(cl.data[7]) / (cl.data[7] + cl.data[4]); + eta.y = static_cast(cl.data[7]) / (cl.data[7] + cl.data[4]); eta.c = cTopRight; break; - // no default to allow compiler to warn about missing cases + // no default to allow compiler to warn about missing cases } return eta; } - -Eta2 calculate_eta2(Cluster2x2 &cl) { +template Eta2 calculate_eta2(Cluster &cl) { Eta2 eta{}; eta.x = static_cast(cl.data[1]) / (cl.data[0] + cl.data[1]); eta.y = static_cast(cl.data[2]) / (cl.data[0] + cl.data[2]); - eta.sum = cl.data[0] + cl.data[1] + cl.data[2]+ cl.data[3]; - eta.c = cBottomLeft; //TODO! This is not correct, but need to put something + eta.sum = cl.data[0] + cl.data[1] + cl.data[2] + cl.data[3]; + eta.c = cBottomLeft; // TODO! This is not correct, but need to put something return eta; } - - int analyze_cluster(Cluster3x3 &cl, int32_t *t2, int32_t *t3, char *quad, double *eta2x, double *eta2y, double *eta3x, double *eta3y) {