cut also for frames

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

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<ROI> m_roi;
-    std::optional<NDArray<int32_t, 2>> m_noise_map;
-    std::optional<NDArray<double, 2>> 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<ROI> m_roi;         /*Region of interest, will be applied if set*/
+    std::optional<NDArray<int32_t, 2>> m_noise_map; /*Noise map to cut photons, will be applied if set*/
+    std::optional<NDArray<double, 2>> m_gain_map; /*Gain map to apply to the clusters, will be applied if set*/
 
   public:
     /**
@@ -98,10 +90,6 @@ class ClusterFile {
      */
     ClusterVector<int32_t> 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<int32_t> &clusters);
     
-    // Need to be migrated to support NDArray and return a ClusterVector
-    // std::vector<Cluster3x3>
-    // 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<int32_t, 2> 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<double, 2> gain_map);
     
     
@@ -137,15 +134,13 @@ class ClusterFile {
     private:
         ClusterVector<int32_t> read_clusters_with_cut(size_t n_clusters);
         ClusterVector<int32_t> read_clusters_without_cut(size_t n_clusters);
+        ClusterVector<int32_t> read_frame_with_cut();
+        ClusterVector<int32_t> 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<double, 2> calculate_eta2(ClusterVector<int> &clusters);
 Eta2 calculate_eta2(Cluster3x3 &cl);
 Eta2 calculate_eta2(Cluster2x2 &cl);

src/ClusterFile.cpp

@@ -60,11 +60,22 @@ void ClusterFile::write_frame(const ClusterVector<int32_t> &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<int32_t> 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<int32_t> ClusterFile::read_frame() {
+ClusterVector<int32_t> 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<int32_t> 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<int32_t> clusters(3, 3, n_clusters);
+    clusters.set_frame_number(frame_number);
+
+    if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) !=
+        static_cast<size_t>(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<int32_t> ClusterFile::read_frame_with_cut() {
     if (m_mode != "r") {
         throw std::runtime_error("File not opened for reading");
     }
@@ -205,22 +259,27 @@ ClusterVector<int32_t> 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<Cluster3x3> clusters(n_clusters);
-    ClusterVector<int32_t> clusters(3, 3, n_clusters);
+    
+    ClusterVector<int32_t> 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<size_t>(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<std::byte*>(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<Cluster3x3> 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<Cluster3x3> 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<void *>(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<void *>(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<double, 2> calculate_eta2(ClusterVector<int> &clusters) {
     //TOTO! make work with 2x2 clusters
     NDArray<double, 2> eta2({static_cast<int64_t>(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<double>(data[4]) / (data[3] + data[4]);
-            if (eta2y && (data[1] + data[4]) != 0)
-                *eta2y = static_cast<double>(data[4]) / (data[1] + data[4]);
-            break;
-        case cBottomRight:
-            if (eta2x && (data[2] + data[5]) != 0)
-                *eta2x = static_cast<double>(data[5]) / (data[4] + data[5]);
-            if (eta2y && (data[1] + data[4]) != 0)
-                *eta2y = static_cast<double>(data[4]) / (data[1] + data[4]);
-            break;
-        case cTopLeft:
-            if (eta2x && (data[7] + data[4]) != 0)
-                *eta2x = static_cast<double>(data[4]) / (data[3] + data[4]);
-            if (eta2y && (data[7] + data[4]) != 0)
-                *eta2y = static_cast<double>(data[7]) / (data[7] + data[4]);
-            break;
-        case cTopRight:
-            if (eta2x && t2max != 0)
-                *eta2x = static_cast<double>(data[5]) / (data[5] + data[4]);
-            if (eta2y && t2max != 0)
-                *eta2y = static_cast<double>(data[7]) / (data[7] + data[4]);
-            break;
-        default:;
-        }
-    }
-
-    if (eta3x || eta3y) {
-        if (eta3x && (data[3] + data[4] + data[5]) != 0)
-            *eta3x = static_cast<double>(-data[3] + data[3 + 2]) /
-                     (data[3] + data[4] + data[5]);
-        if (eta3y && (data[1] + data[4] + data[7]) != 0)
-            *eta3y = static_cast<double>(-data[1] + data[2 * 3 + 1]) /
-                     (data[1] + data[4] + data[7]);
-    }
-
-    return ok;
-}
-
 } // namespace aare

src/ClusterFile.test.cpp

@@ -0,0 +1,65 @@
+#include "aare/ClusterFile.hpp"
+#include "test_config.hpp"
+
+
+#include "aare/defs.hpp"
+#include <catch2/catch_test_macros.hpp>
+#include <filesystem>
+
+
+
+
+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<aare::Cluster3x3>(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);
+}

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");
     }
 }