Interpolate (#137)

- added eta based interpolation

src/ClusterFile.cpp

@@ -108,6 +108,79 @@ ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters) {
     return clusters;
 }
 
+ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters, ROI roi) {
+    if (m_mode != "r") {
+        throw std::runtime_error("File not opened for reading");
+    }
+    
+    ClusterVector<int32_t> clusters(3,3);
+    clusters.reserve(n_clusters);
+
+    int32_t iframe = 0; // frame number needs to be 4 bytes!
+    size_t nph_read = 0;
+    uint32_t nn = m_num_left;
+    uint32_t nph = m_num_left; // number of clusters in frame needs to be 4
+
+    // auto buf = reinterpret_cast<Cluster3x3 *>(clusters.data());
+    // auto buf = clusters.data();
+
+    Cluster3x3 tmp; //this would break if the cluster size changes
+
+    // if there are photons left from previous frame read them first
+    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;
+        }
+        //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<std::byte*>(tmp.data));
+                nph_read++;
+            }
+        }
+
+        m_num_left = nph - nn; // write back the number of photons left
+    }
+
+    if (nph_read < n_clusters) {
+        // keep on reading frames and photons until reaching n_clusters
+        while (fread(&iframe, sizeof(iframe), 1, fp)) {
+            // read number of clusters in frame
+            if (fread(&nph, sizeof(nph), 1, fp)) {
+                if (nph > (n_clusters - nph_read))
+                    nn = n_clusters - nph_read;
+                else
+                    nn = nph;
+
+                // 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<std::byte*>(tmp.data));
+                        nph_read++;
+                    }
+                }
+                m_num_left = nph - nn;
+            }
+            if (nph_read >= n_clusters)
+                break;
+        }
+    }
+
+    // Resize the vector to the number of clusters.
+    // No new allocation, only change bounds.
+    clusters.resize(nph_read);
+    return clusters;
+}
+
 ClusterVector<int32_t> ClusterFile::read_frame() {
     if (m_mode != "r") {
         throw std::runtime_error("File not opened for reading");
@@ -268,11 +341,23 @@ ClusterVector<int32_t> ClusterFile::read_frame() {
 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});
-    for (size_t i = 0; i < clusters.size(); i++) {
-        auto e = calculate_eta2(clusters.at<Cluster3x3>(i));
-        eta2(i, 0) = e.x;
-        eta2(i, 1) = e.y;
+    
+    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<Cluster3x3>(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<Cluster2x2>(i));
+            eta2(i, 0) = e.x;
+            eta2(i, 1) = e.y;
+        }
+    }else{
+        throw std::runtime_error("Only 3x3 and 2x2 clusters are supported");
     }
+    
     return eta2;
 }
 
@@ -290,7 +375,7 @@ Eta2 calculate_eta2(Cluster3x3 &cl) {
     tot2[3] = cl.data[4] + cl.data[5] + cl.data[7] + cl.data[8];
 
     auto c = std::max_element(tot2.begin(), tot2.end()) - tot2.begin();
-
+    eta.sum = tot2[c];
     switch (c) {
     case cBottomLeft:
         if ((cl.data[3] + cl.data[4]) != 0)
@@ -333,6 +418,20 @@ Eta2 calculate_eta2(Cluster3x3 &cl) {
     return eta;
 }
 
+
+Eta2 calculate_eta2(Cluster2x2 &cl) {
+    Eta2 eta{};
+    if ((cl.data[0] + cl.data[1]) != 0)
+        eta.x = static_cast<double>(cl.data[1]) / (cl.data[0] + cl.data[1]);
+    if ((cl.data[0] + cl.data[2]) != 0)
+        eta.y = static_cast<double>(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
+    return eta;
+}
+
+
+
 int analyze_cluster(Cluster3x3 &cl, int32_t *t2, int32_t *t3, char *quad,
                     double *eta2x, double *eta2y, double *eta3x,
                     double *eta3y) {

src/Interpolator.cpp

@@ -0,0 +1,144 @@
+#include "aare/Interpolator.hpp"
+#include "aare/algorithm.hpp"
+
+namespace aare {
+
+Interpolator::Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
+                           NDView<double, 1> ybins, NDView<double, 1> ebins)
+    : m_ietax(etacube), m_ietay(etacube), m_etabinsx(xbins), m_etabinsy(ybins), m_energy_bins(ebins) {
+    if (etacube.shape(0) != xbins.size() || etacube.shape(1) != ybins.size() ||
+        etacube.shape(2) != ebins.size()) {
+        throw std::invalid_argument(
+            "The shape of the etacube does not match the shape of the bins");
+    }
+
+    // Cumulative sum in the x direction
+    for (ssize_t i = 1; i < m_ietax.shape(0); i++) {
+        for (ssize_t j = 0; j < m_ietax.shape(1); j++) {
+            for (ssize_t k = 0; k < m_ietax.shape(2); k++) {
+                m_ietax(i, j, k) += m_ietax(i - 1, j, k);
+            }
+        }
+    }
+
+    // Normalize by the highest row, if norm less than 1 don't do anything
+    for (ssize_t i = 0; i < m_ietax.shape(0); i++) {
+        for (ssize_t j = 0; j < m_ietax.shape(1); j++) {
+            for (ssize_t k = 0; k < m_ietax.shape(2); k++) {
+                auto val = m_ietax(m_ietax.shape(0) - 1, j, k);
+                double norm = val < 1 ? 1 : val;
+                m_ietax(i, j, k) /= norm;
+            }
+        }
+    }
+
+    // Cumulative sum in the y direction
+    for (ssize_t i = 0; i < m_ietay.shape(0); i++) {
+        for (ssize_t j = 1; j < m_ietay.shape(1); j++) {
+            for (ssize_t k = 0; k < m_ietay.shape(2); k++) {
+                m_ietay(i, j, k) += m_ietay(i, j - 1, k);
+            }
+        }
+    }
+
+    // Normalize by the highest column, if norm less than 1 don't do anything
+    for (ssize_t i = 0; i < m_ietay.shape(0); i++) {
+        for (ssize_t j = 0; j < m_ietay.shape(1); j++) {
+            for (ssize_t k = 0; k < m_ietay.shape(2); k++) {
+                auto val = m_ietay(i, m_ietay.shape(1) - 1, k);
+                double norm = val < 1 ? 1 : val;
+                m_ietay(i, j, k) /= norm;
+            }
+        }
+    }
+}
+
+std::vector<Photon> Interpolator::interpolate(const ClusterVector<int32_t>& clusters) {
+    std::vector<Photon> photons;
+    photons.reserve(clusters.size());
+
+    if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) {
+        for (size_t i = 0; i<clusters.size(); i++){
+
+            auto cluster = clusters.at<Cluster3x3>(i);
+            Eta2 eta= calculate_eta2(cluster);
+    
+            Photon photon;
+            photon.x = cluster.x;
+            photon.y = cluster.y;
+            photon.energy = eta.sum;
+    
+            // auto ie = nearest_index(m_energy_bins, photon.energy)-1;
+            // auto ix = nearest_index(m_etabinsx, eta.x)-1;
+            // auto iy = nearest_index(m_etabinsy, eta.y)-1;   
+            //Finding the index of the last element that is smaller
+            //should work fine as long as we have many bins
+            auto ie = last_smaller(m_energy_bins, photon.energy);
+            auto ix = last_smaller(m_etabinsx, eta.x);
+            auto iy = last_smaller(m_etabinsy, eta.y); 
+
+            // fmt::print("ex: {}, ix: {}, iy: {}\n", ie, ix, iy);
+            
+            double dX, dY;
+            int ex, ey;
+            // cBottomLeft = 0,
+            // cBottomRight = 1,
+            // cTopLeft = 2,
+            // cTopRight = 3
+            switch (eta.c) {
+            case cTopLeft:
+                dX = -1.;
+                dY = 0.;
+                break;
+            case cTopRight:;
+                dX = 0.;
+                dY = 0.;
+                break;
+            case cBottomLeft:
+                dX = -1.;
+                dY = -1.;
+                break;
+            case cBottomRight:
+                dX = 0.; 
+                dY = -1.;
+                break;
+            }
+            photon.x += m_ietax(ix, iy, ie)*2 + dX;
+            photon.y += m_ietay(ix, iy, ie)*2 + dY;
+            photons.push_back(photon);
+        }
+    }else if(clusters.cluster_size_x() == 2 || clusters.cluster_size_y() == 2){
+        for (size_t i = 0; i<clusters.size(); i++){
+            auto cluster = clusters.at<Cluster2x2>(i);
+            Eta2 eta= calculate_eta2(cluster);
+    
+            Photon photon;
+            photon.x = cluster.x;
+            photon.y = cluster.y;
+            photon.energy = eta.sum;
+    
+            //Now do some actual interpolation. 
+            //Find which energy bin the cluster is in
+            // auto ie = nearest_index(m_energy_bins, photon.energy)-1;
+            // auto ix = nearest_index(m_etabinsx, eta.x)-1;
+            // auto iy = nearest_index(m_etabinsy, eta.y)-1;  
+            //Finding the index of the last element that is smaller
+            //should work fine as long as we have many bins
+            auto ie = last_smaller(m_energy_bins, photon.energy);
+            auto ix = last_smaller(m_etabinsx, eta.x);
+            auto iy = last_smaller(m_etabinsy, eta.y); 
+
+            photon.x += m_ietax(ix, iy, ie)*2; //eta goes between 0 and 1 but we could move the hit anywhere in the 2x2
+            photon.y += m_ietay(ix, iy, ie)*2;
+            photons.push_back(photon);
+        }
+        
+    }else{
+        throw std::runtime_error("Only 3x3 and 2x2 clusters are supported for interpolation");
+    }
+  
+
+    return photons;
+}
+
+} // namespace aare

src/NDArray.test.cpp

@@ -2,6 +2,7 @@
 #include <array>
 #include <catch2/benchmark/catch_benchmark.hpp>
 #include <catch2/catch_test_macros.hpp>
+#include <numeric>
 
 using aare::NDArray;
 using aare::NDView;
@@ -34,6 +35,24 @@ TEST_CASE("Construct from an NDView") {
     }
 }
 
+TEST_CASE("3D NDArray from NDView"){
+    std::vector<int> data(27);
+    std::iota(data.begin(), data.end(), 0);
+    NDView<int, 3> view(data.data(), Shape<3>{3, 3, 3});
+    NDArray<int, 3> image(view);
+    REQUIRE(image.shape() == view.shape());
+    REQUIRE(image.size() == view.size());
+    REQUIRE(image.data() != view.data());
+
+    for(int64_t i=0; i<image.shape(0); i++){
+        for(int64_t j=0; j<image.shape(1); j++){
+            for(int64_t k=0; k<image.shape(2); k++){
+                REQUIRE(image(i, j, k) == view(i, j, k));
+            }
+        }
+    }
+}
+
 TEST_CASE("1D image") {
     std::array<int64_t, 1> shape{{20}};
     NDArray<short, 1> img(shape, 3);

src/algorithm.test.cpp

@@ -0,0 +1,73 @@
+
+
+#include <catch2/catch_test_macros.hpp>
+#include <aare/algorithm.hpp>
+
+
+TEST_CASE("Find the closed index in a 1D array", "[algorithm]") {
+    aare::NDArray<double, 1> arr({5});
+    for (size_t i = 0; i < arr.size(); i++) {
+        arr[i] = i;
+    }
+    // arr 0, 1, 2, 3, 4
+    REQUIRE(aare::nearest_index(arr, 2.3) == 2);
+    REQUIRE(aare::nearest_index(arr, 2.6) == 3);
+    REQUIRE(aare::nearest_index(arr, 45.0) == 4);
+    REQUIRE(aare::nearest_index(arr, 0.0) == 0);
+    REQUIRE(aare::nearest_index(arr, -1.0) == 0);
+}
+
+TEST_CASE("Passing integers to nearest_index works", "[algorithm]"){
+    aare::NDArray<int, 1> arr({5});
+    for (size_t i = 0; i < arr.size(); i++) {
+        arr[i] = i;
+    }
+    // arr 0, 1, 2, 3, 4
+    REQUIRE(aare::nearest_index(arr, 2) == 2);
+    REQUIRE(aare::nearest_index(arr, 3) == 3);
+    REQUIRE(aare::nearest_index(arr, 45) == 4);
+    REQUIRE(aare::nearest_index(arr, 0) == 0);
+    REQUIRE(aare::nearest_index(arr, -1) == 0);
+}
+
+
+TEST_CASE("nearest_index works with std::vector", "[algorithm]"){
+    std::vector<double> vec = {0, 1, 2, 3, 4};
+    REQUIRE(aare::nearest_index(vec, 2.123) == 2);
+    REQUIRE(aare::nearest_index(vec, 2.66) == 3);
+    REQUIRE(aare::nearest_index(vec, 4555555.0) == 4);
+    REQUIRE(aare::nearest_index(vec, 0.0) == 0);
+    REQUIRE(aare::nearest_index(vec, -10.0) == 0);
+}
+
+TEST_CASE("nearest index works with std::array", "[algorithm]"){
+    std::array<double, 5> arr = {0, 1, 2, 3, 4};
+    REQUIRE(aare::nearest_index(arr, 2.123) == 2);
+    REQUIRE(aare::nearest_index(arr, 2.501) == 3);
+    REQUIRE(aare::nearest_index(arr, 4555555.0) == 4);
+    REQUIRE(aare::nearest_index(arr, 0.0) == 0);
+    REQUIRE(aare::nearest_index(arr, -10.0) == 0);
+}
+
+
+TEST_CASE("last smaller", "[algorithm]"){
+    aare::NDArray<double, 1> arr({5});
+    for (size_t i = 0; i < arr.size(); i++) {
+        arr[i] = i;
+    }
+    // arr 0, 1, 2, 3, 4
+    REQUIRE(aare::last_smaller(arr, -10.0) == 0);
+    REQUIRE(aare::last_smaller(arr, 0.0) == 0);
+    REQUIRE(aare::last_smaller(arr, 2.3) == 2);
+    REQUIRE(aare::last_smaller(arr, 253.) == 4);
+}
+
+TEST_CASE("returns last bin strictly smaller", "[algorithm]"){
+    aare::NDArray<double, 1> arr({5});
+    for (size_t i = 0; i < arr.size(); i++) {
+        arr[i] = i;
+    }
+    // arr 0, 1, 2, 3, 4
+    REQUIRE(aare::last_smaller(arr, 2.0) == 2);
+
+}