diff --git a/CMakeLists.txt b/CMakeLists.txt index 4772f0b..6275e13 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -396,6 +396,7 @@ if(AARE_TESTS) ${CMAKE_CURRENT_SOURCE_DIR}/src/NDView.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinder.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterVector.test.cpp + ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFile.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/Pedestal.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.test.cpp ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.test.cpp diff --git a/include/aare/ClusterFile.hpp b/include/aare/ClusterFile.hpp index 87f5783..22f4183 100644 --- a/include/aare/ClusterFile.hpp +++ b/include/aare/ClusterFile.hpp @@ -10,7 +10,7 @@ namespace aare { - +//TODO! Legacy enums, migrate to enum class typedef enum { cBottomLeft = 0, cBottomRight = 1, @@ -44,15 +44,7 @@ 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 -.... -*/ + /** * @brief Class to read and write cluster files @@ -61,19 +53,19 @@ uint32_t number_of_clusters * * int32_t frame_number * uint32_t number_of_clusters - * int16_t x, int16_t y, int32_t data[9] x number_of_clusters + * int16_t x, int16_t y, int32_t data[9] * number_of_clusters * int32_t frame_number * uint32_t number_of_clusters * etc. */ class ClusterFile { FILE *fp{}; - uint32_t m_num_left{}; - size_t m_chunk_size{}; - const std::string m_mode; - std::optional m_roi; - std::optional> m_noise_map; - std::optional> m_gain_map; + uint32_t m_num_left{}; /*Number of photons left in frame*/ + size_t m_chunk_size{}; /*Number of clusters to read at a time*/ + const std::string m_mode; /*Mode to open the file in*/ + std::optional m_roi; /*Region of interest, will be applied if set*/ + std::optional> m_noise_map; /*Noise map to cut photons, will be applied if set*/ + std::optional> m_gain_map; /*Gain map to apply to the clusters, will be applied if set*/ public: /** @@ -98,10 +90,6 @@ class ClusterFile { */ ClusterVector read_clusters(size_t n_clusters); - - - - /** * @brief Read a single frame from the file and return the clusters. The * cluster vector will have the frame number set. @@ -113,19 +101,28 @@ class ClusterFile { 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); - /** * @brief Return the chunk size */ size_t chunk_size() const { return m_chunk_size; } + /** + * @brief Set the region of interest to use when reading clusters. If set only clusters within + * the ROI will be read. + */ void set_roi(ROI roi); + /** + * @brief Set the noise map to use when reading clusters. If set clusters below the noise + * level will be discarded. Selection criteria one of: Central pixel above noise, highest + * 2x2 sum above 2 * noise, total sum above 3 * noise. + */ void set_noise_map(const NDView noise_map); + /** + * @brief Set the gain map to use when reading clusters. If set the gain map will be applied + * to the clusters that pass ROI and noise_map selection. + */ void set_gain_map(const NDView gain_map); @@ -137,15 +134,13 @@ class ClusterFile { private: ClusterVector read_clusters_with_cut(size_t n_clusters); ClusterVector read_clusters_without_cut(size_t n_clusters); + ClusterVector read_frame_with_cut(); + ClusterVector read_frame_without_cut(); bool is_selected(Cluster3x3 &cl); Cluster3x3 read_one_cluster(); }; -int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, - double *eta2x, double *eta2y, double *eta3x, double *eta3y); -int analyze_cluster(Cluster3x3 &cl, int32_t *t2, int32_t *t3, char *quad, - double *eta2x, double *eta2y, double *eta3x, double *eta3y); - +//TODO! helper functions that doesn't really belong here NDArray calculate_eta2(ClusterVector &clusters); Eta2 calculate_eta2(Cluster3x3 &cl); Eta2 calculate_eta2(Cluster2x2 &cl); diff --git a/src/ClusterFile.cpp b/src/ClusterFile.cpp index 99e2f60..f4ef0ae 100644 --- a/src/ClusterFile.cpp +++ b/src/ClusterFile.cpp @@ -60,11 +60,22 @@ void ClusterFile::write_frame(const ClusterVector &clusters) { !(clusters.cluster_size_y() == 3)) { throw std::runtime_error("Only 3x3 clusters are supported"); } + //First write the frame number - 4 bytes int32_t frame_number = clusters.frame_number(); - fwrite(&frame_number, sizeof(frame_number), 1, fp); + if(fwrite(&frame_number, sizeof(frame_number), 1, fp)!=1){ + throw std::runtime_error(LOCATION + "Could not write frame number"); + } + + //Then write the number of clusters - 4 bytes uint32_t n_clusters = clusters.size(); - fwrite(&n_clusters, sizeof(n_clusters), 1, fp); - fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp); + if(fwrite(&n_clusters, sizeof(n_clusters), 1, fp)!=1){ + throw std::runtime_error(LOCATION + "Could not write number of clusters"); + } + + //Now write the clusters in the frame + if(fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp)!=clusters.size()){ + throw std::runtime_error(LOCATION + "Could not write clusters"); + } } @@ -110,6 +121,7 @@ ClusterVector ClusterFile::read_clusters_without_cut(size_t n_clusters) if (nph_read < n_clusters) { // keep on reading frames and photons until reaching n_clusters while (fread(&iframe, sizeof(iframe), 1, fp)) { + clusters.set_frame_number(iframe); // read number of clusters in frame if (fread(&nph, sizeof(nph), 1, fp)) { if (nph > (n_clusters - nph_read)) @@ -192,7 +204,49 @@ Cluster3x3 ClusterFile::read_one_cluster(){ return c; } -ClusterVector ClusterFile::read_frame() { +ClusterVector ClusterFile::read_frame(){ + if (m_mode != "r") { + throw std::runtime_error(LOCATION + "File not opened for reading"); + } + if (m_noise_map || m_roi){ + return read_frame_with_cut(); + }else{ + return read_frame_without_cut(); + } +} + +ClusterVector ClusterFile::read_frame_without_cut() { + 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"); + } + int32_t frame_number; + if (fread(&frame_number, sizeof(frame_number), 1, fp) != 1) { + throw std::runtime_error(LOCATION + "Could not read frame number"); + } + + int32_t n_clusters; // Saved as 32bit integer in the cluster file + if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) { + throw std::runtime_error(LOCATION + "Could not read number of clusters"); + } + + ClusterVector clusters(3, 3, n_clusters); + clusters.set_frame_number(frame_number); + + if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) != + static_cast(n_clusters)) { + throw std::runtime_error(LOCATION + "Could not read clusters"); + } + clusters.resize(n_clusters); + if (m_gain_map) + clusters.apply_gain_map(m_gain_map->view()); + return clusters; +} + +ClusterVector ClusterFile::read_frame_with_cut() { if (m_mode != "r") { throw std::runtime_error("File not opened for reading"); } @@ -205,22 +259,27 @@ ClusterVector ClusterFile::read_frame() { throw std::runtime_error("Could not read frame number"); } - int32_t n_clusters; // Saved as 32bit integer in the cluster file - if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) { + + if (fread(&m_num_left, sizeof(m_num_left), 1, fp) != 1) { throw std::runtime_error("Could not read number of clusters"); } - // std::vector clusters(n_clusters); - ClusterVector clusters(3, 3, n_clusters); + + ClusterVector clusters(3, 3); + clusters.reserve(m_num_left); clusters.set_frame_number(frame_number); - - if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) != - static_cast(n_clusters)) { - throw std::runtime_error("Could not read clusters"); + while(m_num_left){ + Cluster3x3 c = read_one_cluster(); + if(is_selected(c)){ + clusters.push_back(c.x, c.y, reinterpret_cast(c.data)); + } } - clusters.resize(n_clusters); + if (m_gain_map) + clusters.apply_gain_map(m_gain_map->view()); return clusters; } + + bool ClusterFile::is_selected(Cluster3x3 &cl) { //Should fail fast if (m_roi) { @@ -242,133 +301,6 @@ bool ClusterFile::is_selected(Cluster3x3 &cl) { return true; } -// 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 -// // nx, int ny) { -// int iframe = 0; -// // uint32_t nph = *n_left; -// uint32_t nph = m_num_left; -// // uint32_t nn = *n_left; -// uint32_t nn = m_num_left; -// size_t nph_read = 0; - -// int32_t t2max, tot1; -// int32_t tot3; -// // Cluster *ptr = buf; -// Cluster3x3 *ptr = clusters.data(); -// int good = 1; -// double noise; -// // read photons left from previous frame -// if (noise_map) -// printf("Using noise map\n"); - -// if (nph) { -// if (nph > n_clusters) { -// // if we have more photons left in the frame then photons to -// // read we read directly the requested number -// nn = n_clusters; -// } else { -// nn = nph; -// } -// 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); -// if (n_read != 1) { -// clusters.resize(nph_read); -// return clusters; -// } -// // TODO! error handling on 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]; -// if (tot1 > noise || t2max > 2 * noise || tot3 > 3 * noise) { -// ; -// } else { -// good = 0; -// printf("%d %d %f %d %d %d\n", ptr->x, ptr->y, noise, -// tot1, t2max, tot3); -// } -// } 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) { -// // // 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(Cluster3x3), 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; -// } -// } -// // printf("%d\n",nph_read); -// clusters.resize(nph_read); -// return clusters; -// } - NDArray calculate_eta2(ClusterVector &clusters) { //TOTO! make work with 2x2 clusters NDArray eta2({static_cast(clusters.size()), 2}); @@ -462,111 +394,4 @@ Eta2 calculate_eta2(Cluster2x2 &cl) { } - -int analyze_cluster(Cluster3x3 &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 analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad, - double *eta2x, double *eta2y, double *eta3x, double *eta3y) { - - int ok = 1; - - int32_t tot2[4]; - int32_t t2max = 0; - char c = 0; - int32_t val, tot3; - - tot3 = 0; - for (int i = 0; i < 4; i++) - tot2[i] = 0; - - for (int ix = 0; ix < 3; ix++) { - for (int iy = 0; iy < 3; iy++) { - val = data[iy * 3 + ix]; - // printf ("%d ",data[iy * 3 + ix]); - tot3 += val; - if (ix <= 1 && iy <= 1) - tot2[cBottomLeft] += val; - if (ix >= 1 && iy <= 1) - tot2[cBottomRight] += val; - if (ix <= 1 && iy >= 1) - tot2[cTopLeft] += val; - if (ix >= 1 && iy >= 1) - tot2[cTopRight] += val; - } - // printf ("\n"); - } - // printf ("\n"); - - if (t2 || quad) { - - t2max = tot2[0]; - c = cBottomLeft; - for (int i = 1; i < 4; i++) { - if (tot2[i] > t2max) { - t2max = tot2[i]; - c = i; - } - } - // printf("*** %d %d %d %d -- - // %d\n",tot2[0],tot2[1],tot2[2],tot2[3],t2max); - if (quad) - *quad = c; - if (t2) - *t2 = t2max; - } - - if (t3) - *t3 = tot3; - - if (eta2x || eta2y) { - if (eta2x) - *eta2x = 0; - if (eta2y) - *eta2y = 0; - switch (c) { - case cBottomLeft: - if (eta2x && (data[3] + data[4]) != 0) - *eta2x = static_cast(data[4]) / (data[3] + data[4]); - if (eta2y && (data[1] + data[4]) != 0) - *eta2y = static_cast(data[4]) / (data[1] + data[4]); - break; - case cBottomRight: - if (eta2x && (data[2] + data[5]) != 0) - *eta2x = static_cast(data[5]) / (data[4] + data[5]); - if (eta2y && (data[1] + data[4]) != 0) - *eta2y = static_cast(data[4]) / (data[1] + data[4]); - break; - case cTopLeft: - if (eta2x && (data[7] + data[4]) != 0) - *eta2x = static_cast(data[4]) / (data[3] + data[4]); - if (eta2y && (data[7] + data[4]) != 0) - *eta2y = static_cast(data[7]) / (data[7] + data[4]); - break; - case cTopRight: - if (eta2x && t2max != 0) - *eta2x = static_cast(data[5]) / (data[5] + data[4]); - if (eta2y && t2max != 0) - *eta2y = static_cast(data[7]) / (data[7] + data[4]); - break; - default:; - } - } - - if (eta3x || eta3y) { - if (eta3x && (data[3] + data[4] + data[5]) != 0) - *eta3x = static_cast(-data[3] + data[3 + 2]) / - (data[3] + data[4] + data[5]); - if (eta3y && (data[1] + data[4] + data[7]) != 0) - *eta3y = static_cast(-data[1] + data[2 * 3 + 1]) / - (data[1] + data[4] + data[7]); - } - - return ok; -} - } // namespace aare \ No newline at end of file diff --git a/src/ClusterFile.test.cpp b/src/ClusterFile.test.cpp new file mode 100644 index 0000000..8b2dc8d --- /dev/null +++ b/src/ClusterFile.test.cpp @@ -0,0 +1,65 @@ +#include "aare/ClusterFile.hpp" +#include "test_config.hpp" + + +#include "aare/defs.hpp" +#include +#include + + + + +using aare::ClusterFile; + +TEST_CASE("Read one frame from a a cluster file", "[.integration]") { + //We know that the frame has 97 clusters + auto fpath = test_data_path() / "clusters" / "single_frame_97_clustrers.clust"; + REQUIRE(std::filesystem::exists(fpath)); + + ClusterFile f(fpath); + auto clusters = f.read_frame(); + REQUIRE(clusters.size() == 97); + REQUIRE(clusters.frame_number() == 135); +} + +TEST_CASE("Read one frame using ROI", "[.integration]") { + //We know that the frame has 97 clusters + auto fpath = test_data_path() / "clusters" / "single_frame_97_clustrers.clust"; + REQUIRE(std::filesystem::exists(fpath)); + + ClusterFile f(fpath); + aare::ROI roi; + roi.xmin = 0; + roi.xmax = 50; + roi.ymin = 200; + roi.ymax = 249; + f.set_roi(roi); + auto clusters = f.read_frame(); + REQUIRE(clusters.size() == 49); + REQUIRE(clusters.frame_number() == 135); + + //Check that all clusters are within the ROI + for (size_t i = 0; i < clusters.size(); i++) { + auto c = clusters.at(i); + REQUIRE(c.x >= roi.xmin); + REQUIRE(c.x <= roi.xmax); + REQUIRE(c.y >= roi.ymin); + REQUIRE(c.y <= roi.ymax); + } + +} + + +TEST_CASE("Read clusters from single frame file", "[.integration]") { + + auto fpath = test_data_path() / "clusters" / "single_frame_97_clustrers.clust"; + REQUIRE(std::filesystem::exists(fpath)); + + ClusterFile f(fpath); + auto clusters = f.read_clusters(500); + REQUIRE(clusters.size() == 97); + + + //Cluster vector should hold the last read frame number: + REQUIRE(clusters.frame_number() == 135); +} diff --git a/tests/test_config.hpp.in b/tests/test_config.hpp.in index 62993b7..e314b8f 100644 --- a/tests/test_config.hpp.in +++ b/tests/test_config.hpp.in @@ -7,6 +7,6 @@ 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"); + throw std::runtime_error("Path to test data: $AARE_TEST_DATA not set"); } } \ No newline at end of file