templated calculate_pedestal with boolean template argument only_gain… (#218)

some refactoring for less code duplication, added functionality
drop_dimension in NDArray

CMakeLists.txt

@@ -438,6 +438,7 @@ endif()
 if(AARE_TESTS)
     set(TestSources
             ${CMAKE_CURRENT_SOURCE_DIR}/src/algorithm.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/calibration.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/defs.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/decode.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.test.cpp

include/aare/NDArray.hpp

@@ -25,7 +25,7 @@ template <typename T, ssize_t Ndim = 2>
 class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
     std::array<ssize_t, Ndim> shape_;
     std::array<ssize_t, Ndim> strides_;
-    size_t size_{};
+    size_t size_{}; //TODO! do we need to store size when we have shape?
     T *data_;
 
   public:
@@ -33,7 +33,7 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
      * @brief Default constructor. Will construct an empty NDArray.
      *
      */
-    NDArray() : shape_(), strides_(c_strides<Ndim>(shape_)), data_(nullptr){};
+    NDArray() : shape_(), strides_(c_strides<Ndim>(shape_)), data_(nullptr) {};
 
     /**
      * @brief Construct a new NDArray object with a given shape.
@@ -43,8 +43,7 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
      */
     explicit NDArray(std::array<ssize_t, Ndim> shape)
         : shape_(shape), strides_(c_strides<Ndim>(shape_)),
-          size_(std::accumulate(shape_.begin(), shape_.end(), 1,
+          size_(num_elements(shape_)),
-                                std::multiplies<>())),
           data_(new T[size_]) {}
 
     /**
@@ -79,6 +78,24 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
         other.reset(); // TODO! is this necessary?
     }
 
+
+    //Move constructor from an an array with Ndim + 1
+    template <ssize_t M, typename = std::enable_if_t<(M == Ndim + 1)>>
+    NDArray(NDArray<T, M> &&other) 
+        : shape_(drop_first_dim(other.shape())),
+          strides_(c_strides<Ndim>(shape_)), size_(num_elements(shape_)),
+          data_(other.data()) {
+
+            // For now only allow move if the size matches, to avoid unreachable data
+            // if the use case arises we can remove this check
+            if(size() != other.size()) {
+                data_ = nullptr; // avoid double free, other will clean up the memory in it's destructor
+                throw std::runtime_error(LOCATION +
+                                         "Size mismatch in move constructor of NDArray<T, Ndim-1>");
+            }
+        other.reset();
+    }
+
     // Copy constructor
     NDArray(const NDArray &other)
         : shape_(other.shape_), strides_(c_strides<Ndim>(shape_)),
@@ -217,7 +234,6 @@ class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
         std::fill(shape_.begin(), shape_.end(), 0);
         std::fill(strides_.begin(), strides_.end(), 0);
     }
-
 };
 
 // Move assign
@@ -382,8 +398,6 @@ NDArray<T, Ndim> NDArray<T, Ndim>::operator*(const T &value) {
     return result;
 }
 
-
-
 template <typename T, ssize_t Ndim>
 std::ostream &operator<<(std::ostream &os, const NDArray<T, Ndim> &arr) {
     for (auto row = 0; row < arr.shape(0); ++row) {
@@ -434,17 +448,18 @@ NDArray<T, Ndim> load(const std::string &pathname,
     return img;
 }
 
-
 template <typename RT, typename NT, typename DT, ssize_t Ndim>
 NDArray<RT, Ndim> safe_divide(const NDArray<NT, Ndim> &numerator,
-                                 const NDArray<DT, Ndim> &denominator) {
+                              const NDArray<DT, Ndim> &denominator) {
     if (numerator.shape() != denominator.shape()) {
-        throw std::runtime_error("Shapes of numerator and denominator must match");
+        throw std::runtime_error(
+            "Shapes of numerator and denominator must match");
     }
     NDArray<RT, Ndim> result(numerator.shape());
     for (ssize_t i = 0; i < numerator.size(); ++i) {
         if (denominator[i] != 0) {
-            result[i] = static_cast<RT>(numerator[i]) / static_cast<RT>(denominator[i]);
+            result[i] =
+                static_cast<RT>(numerator[i]) / static_cast<RT>(denominator[i]);
         } else {
             result[i] = RT{0}; // or handle division by zero as needed
         }

include/aare/NDView.hpp

@@ -26,6 +26,33 @@ Shape<Ndim> make_shape(const std::vector<size_t> &shape) {
     return arr;
 }
 
+
+/**
+ * @brief Helper function to drop the first dimension of a shape.
+ * This is useful when you want to create a 2D view from a 3D array.
+ * @param shape The shape to drop the first dimension from.
+ * @return A new shape with the first dimension dropped.
+ */
+template<size_t Ndim>
+Shape<Ndim-1> drop_first_dim(const Shape<Ndim> &shape) {
+    static_assert(Ndim > 1, "Cannot drop first dimension from a 1D shape");
+    Shape<Ndim - 1> new_shape;
+    std::copy(shape.begin() + 1, shape.end(), new_shape.begin());
+    return new_shape;
+}
+
+/**
+ * @brief Helper function when constructing NDArray/NDView. Calculates the number
+ * of elements in the resulting array from a shape.
+ * @param shape The shape to calculate the number of elements for.
+ * @return The number of elements in and NDArray/NDView of that shape.
+ */
+template <size_t Ndim>
+size_t num_elements(const Shape<Ndim> &shape) {
+    return std::accumulate(shape.begin(), shape.end(), 1,
+                           std::multiplies<size_t>());
+}
+
 template <ssize_t Dim = 0, typename Strides>
 ssize_t element_offset(const Strides & /*unused*/) {
     return 0;

include/aare/calibration.hpp

@@ -103,7 +103,7 @@ void apply_calibration_impl(NDView<T, 3> res, NDView<uint16_t, 3> raw_data,
     }
 }
 
-template <class T, ssize_t Ndim=3>
+template <class T, ssize_t Ndim = 3>
 void apply_calibration(NDView<T, 3> res, NDView<uint16_t, 3> raw_data,
                        NDView<T, Ndim> ped, NDView<T, Ndim> cal,
                        ssize_t n_threads = 4) {
@@ -113,52 +113,44 @@ void apply_calibration(NDView<T, 3> res, NDView<uint16_t, 3> raw_data,
     for (const auto &lim : limits)
         futures.push_back(std::async(
             static_cast<void (*)(NDView<T, 3>, NDView<uint16_t, 3>,
-                                 NDView<T, Ndim>, NDView<T, Ndim>, int,
+                                 NDView<T, Ndim>, NDView<T, Ndim>, int, int)>(
-                                 int)>(
                 apply_calibration_impl),
             res, raw_data, ped, cal, lim.first, lim.second));
     for (auto &f : futures)
         f.get();
 }
 
-
+template <bool only_gain0>
 std::pair<NDArray<size_t, 3>, NDArray<size_t, 3>>
-sum_and_count_per_gain(NDView<uint16_t, 3> raw_data);
+sum_and_count_per_gain(NDView<uint16_t, 3> raw_data) {
-
+    constexpr ssize_t num_gains = only_gain0 ? 1 : 3;
-std::pair<NDArray<size_t, 2>, NDArray<size_t, 2>>
+    NDArray<size_t, 3> accumulator(
-sum_and_count_g0(NDView<uint16_t, 3> raw_data);
+        std::array<ssize_t, 3>{num_gains, raw_data.shape(1), raw_data.shape(2)},
-
+        0);
-template <typename T>
+    NDArray<size_t, 3> count(
-NDArray<T, 2> calculate_pedestal_g0(NDView<uint16_t, 3> raw_data,
+        std::array<ssize_t, 3>{num_gains, raw_data.shape(1), raw_data.shape(2)},
-                                    ssize_t n_threads) {
+        0);
-    std::vector<std::future<std::pair<NDArray<size_t, 2>, NDArray<size_t, 2>>>>
+    for (int frame_nr = 0; frame_nr != raw_data.shape(0); ++frame_nr) {
-        futures;
+        for (int row = 0; row != raw_data.shape(1); ++row) {
-    futures.reserve(n_threads);
+            for (int col = 0; col != raw_data.shape(2); ++col) {
-
+                auto [value, gain] =
-    auto subviews = make_subviews(raw_data, n_threads);
+                    get_value_and_gain(raw_data(frame_nr, row, col));
-
+                if (gain != 0 && only_gain0)
-    for (auto view : subviews) {
+                    continue;
-        futures.push_back(std::async(
+                accumulator(gain, row, col) += value;
-            static_cast<std::pair<NDArray<size_t, 2>, NDArray<size_t, 2>> (*)(
+                count(gain, row, col) += 1;
-                NDView<uint16_t, 3>)>(&sum_and_count_g0),
+            }
-            view));
+        }
-    }
-    NDArray<size_t, 2> accumulator(
-        std::array<ssize_t, 2>{raw_data.shape(1), raw_data.shape(2)}, 0);
-    NDArray<size_t, 2> count(
-        std::array<ssize_t, 2>{raw_data.shape(1), raw_data.shape(2)}, 0);
-    for (auto &f : futures) {
-        auto [acc, cnt] = f.get();
-        accumulator += acc;
-        count += cnt;
     }
 
-    return safe_divide<T>(accumulator, count);
+    return {std::move(accumulator), std::move(count)};
 }
 
-template <typename T>
+template <typename T, bool only_gain0 = false>
-NDArray<T, 3> calculate_pedestal(NDView<uint16_t, 3> raw_data,
+NDArray<T, 3 - static_cast<ssize_t>(only_gain0)>
-                                 ssize_t n_threads) {
+calculate_pedestal(NDView<uint16_t, 3> raw_data, ssize_t n_threads) {
+
+    constexpr ssize_t num_gains = only_gain0 ? 1 : 3;
     std::vector<std::future<std::pair<NDArray<size_t, 3>, NDArray<size_t, 3>>>>
         futures;
     futures.reserve(n_threads);
@@ -168,21 +160,25 @@ NDArray<T, 3> calculate_pedestal(NDView<uint16_t, 3> raw_data,
     for (auto view : subviews) {
         futures.push_back(std::async(
             static_cast<std::pair<NDArray<size_t, 3>, NDArray<size_t, 3>> (*)(
-                NDView<uint16_t, 3>)>(&sum_and_count_per_gain),
+                NDView<uint16_t, 3>)>(&sum_and_count_per_gain<only_gain0>),
             view));
     }
+    Shape<3> shape{num_gains, raw_data.shape(1), raw_data.shape(2)};
+    NDArray<size_t, 3> accumulator(shape, 0);
+    NDArray<size_t, 3> count(shape, 0);
 
-    NDArray<size_t, 3> accumulator(
+    // Combine the results from the futures
-        std::array<ssize_t, 3>{3, raw_data.shape(1), raw_data.shape(2)}, 0);
-    NDArray<size_t, 3> count(
-        std::array<ssize_t, 3>{3, raw_data.shape(1), raw_data.shape(2)}, 0);
     for (auto &f : futures) {
         auto [acc, cnt] = f.get();
         accumulator += acc;
         count += cnt;
     }
 
+
+    // Will move to a NDArray<T, 3 - static_cast<ssize_t>(only_gain0)>
+    // if only_gain0 is true
     return safe_divide<T>(accumulator, count);
+
 }
 
 /**
@@ -205,4 +201,9 @@ NDArray<int, 2> count_switching_pixels(NDView<uint16_t, 3> raw_data);
 NDArray<int, 2> count_switching_pixels(NDView<uint16_t, 3> raw_data,
                                        ssize_t n_threads);
 
+template <typename T>
+auto calculate_pedestal_g0(NDView<uint16_t, 3> raw_data, ssize_t n_threads) {
+    return calculate_pedestal<T, true>(raw_data, n_threads);
+}
+
 } // namespace aare

python/src/bind_calibration.hpp

@@ -56,7 +56,7 @@ py::array_t<T> pybind_calculate_pedestal(
 
     auto data_span = make_view_3d(data);
     auto arr = new NDArray<T, 3>{};
-    *arr = aare::calculate_pedestal<T>(data_span, n_threads);
+    *arr = aare::calculate_pedestal<T, false>(data_span, n_threads);
     return return_image_data(arr);
 }
 
@@ -67,7 +67,7 @@ py::array_t<T> pybind_calculate_pedestal_g0(
 
     auto data_span = make_view_3d(data);
     auto arr = new NDArray<T, 2>{};
-    *arr = aare::calculate_pedestal_g0<T>(data_span, n_threads);
+    *arr = aare::calculate_pedestal<T, true>(data_span, n_threads);
     return return_image_data(arr);
 }
 

src/NDArray.test.cpp

@@ -428,3 +428,29 @@ TEST_CASE("Construct an NDArray from an std::array") {
         REQUIRE(a(i) == b[i]);
     }
 }
+
+
+
+TEST_CASE("Move construct from an array with Ndim + 1") {
+    NDArray<int, 3> a({{1,2,2}}, 0);
+    a(0, 0, 0) = 1;
+    a(0, 0, 1) = 2;
+    a(0, 1, 0) = 3;
+    a(0, 1, 1) = 4;
+
+
+    NDArray<int, 2> b(std::move(a));
+    REQUIRE(b.shape() == Shape<2>{2,2});
+    REQUIRE(b.size() == 4);
+    REQUIRE(b(0, 0) == 1);
+    REQUIRE(b(0, 1) == 2);
+    REQUIRE(b(1, 0) == 3);
+    REQUIRE(b(1, 1) == 4);
+
+}
+
+TEST_CASE("Move construct from an array with Ndim + 1 throws on size mismatch") {
+    NDArray<int, 3> a({{2,2,2}}, 0);
+    REQUIRE_THROWS(NDArray<int, 2>(std::move(a)));
+}
+

src/calibration.cpp

@@ -2,46 +2,14 @@
 
 namespace aare {
 
-std::pair<NDArray<size_t, 3>, NDArray<size_t,3>> sum_and_count_per_gain(NDView<uint16_t,3> raw_data){
-    NDArray<size_t, 3> accumulator(std::array<ssize_t, 3>{3, raw_data.shape(1), raw_data.shape(2)}, 0);
-    NDArray<size_t, 3> count(std::array<ssize_t, 3>{3, raw_data.shape(1), raw_data.shape(2)}, 0);
-    for (int frame_nr = 0; frame_nr != raw_data.shape(0); ++frame_nr) {
-        for (int row = 0; row != raw_data.shape(1); ++row) {
-            for (int col = 0; col != raw_data.shape(2); ++col) {
-                auto [value, gain] = get_value_and_gain(raw_data(frame_nr, row, col));
-                accumulator(gain, row, col) += value;
-                count(gain, row, col) += 1;
-            }
-        }
-    }
-
-    return {std::move(accumulator), std::move(count)};
-}
-
-std::pair<NDArray<size_t, 2>, NDArray<size_t,2>> sum_and_count_g0(NDView<uint16_t, 3> raw_data){
-    NDArray<size_t, 2> accumulator(std::array<ssize_t, 2>{raw_data.shape(1), raw_data.shape(2)}, 0);
-    NDArray<size_t, 2> count(std::array<ssize_t, 2>{raw_data.shape(1), raw_data.shape(2)}, 0);
-    for (int frame_nr = 0; frame_nr != raw_data.shape(0); ++frame_nr) {
-        for (int row = 0; row != raw_data.shape(1); ++row) {
-            for (int col = 0; col != raw_data.shape(2); ++col) {
-                auto [value, gain] = get_value_and_gain(raw_data(frame_nr, row, col));
-                if (gain != 0)
-                    continue; // we only care about gain 0
-                accumulator(row, col) += value;
-                count(row, col) += 1;
-            }
-        }
-    }
-
-    return {std::move(accumulator), std::move(count)};
-}
-
 NDArray<int, 2> count_switching_pixels(NDView<uint16_t, 3> raw_data) {
-    NDArray<int, 2> switched(std::array<ssize_t, 2>{raw_data.shape(1), raw_data.shape(2)},0);
+    NDArray<int, 2> switched(
+        std::array<ssize_t, 2>{raw_data.shape(1), raw_data.shape(2)}, 0);
     for (int frame_nr = 0; frame_nr != raw_data.shape(0); ++frame_nr) {
         for (int row = 0; row != raw_data.shape(1); ++row) {
             for (int col = 0; col != raw_data.shape(2); ++col) {
-                auto [value, gain] = get_value_and_gain(raw_data(frame_nr, row, col));
+                auto [value, gain] =
+                    get_value_and_gain(raw_data(frame_nr, row, col));
                 if (gain != 0) {
                     switched(row, col) += 1;
                 }
@@ -51,16 +19,20 @@ NDArray<int, 2> count_switching_pixels(NDView<uint16_t, 3> raw_data) {
     return switched;
 }
 
-
+NDArray<int, 2> count_switching_pixels(NDView<uint16_t, 3> raw_data,
-NDArray<int, 2> count_switching_pixels(NDView<uint16_t, 3> raw_data, ssize_t n_threads){
+                                       ssize_t n_threads) {
-    NDArray<int, 2> switched(std::array<ssize_t, 2>{raw_data.shape(1), raw_data.shape(2)},0);
+    NDArray<int, 2> switched(
+        std::array<ssize_t, 2>{raw_data.shape(1), raw_data.shape(2)}, 0);
     std::vector<std::future<NDArray<int, 2>>> futures;
     futures.reserve(n_threads);
 
     auto subviews = make_subviews(raw_data, n_threads);
 
     for (auto view : subviews) {
-        futures.push_back(std::async(static_cast<NDArray<int, 2>(*)(NDView<uint16_t, 3>)>(&count_switching_pixels), view));
+        futures.push_back(
+            std::async(static_cast<NDArray<int, 2> (*)(NDView<uint16_t, 3>)>(
+                           &count_switching_pixels),
+                       view));
     }
 
     for (auto &f : futures) {
@@ -69,4 +41,4 @@ NDArray<int, 2> count_switching_pixels(NDView<uint16_t, 3> raw_data, ssize_t n_t
     return switched;
 }
 
-}// namespace aare
+} // namespace aare

src/calibration.test.cpp

@@ -0,0 +1,49 @@
+/************************************************
+ * @file test-Cluster.cpp
+ * @short test case for generic Cluster, ClusterVector, and calculate_eta2
+ ***********************************************/
+
+#include "aare/calibration.hpp"
+
+// #include "catch.hpp"
+#include <array>
+#include <catch2/catch_all.hpp>
+#include <catch2/catch_test_macros.hpp>
+
+using namespace aare;
+
+TEST_CASE("Test Pedestal Generation", "[.calibration]") {
+    NDArray<uint16_t, 3> raw(std::array<ssize_t, 3>{3, 2, 2}, 0);
+
+    // gain 0
+    raw(0, 0, 0) = 100;
+    raw(1, 0, 0) = 200;
+    raw(2, 0, 0) = 300;
+
+    // gain 1
+    raw(0, 0, 1) = (1 << 14) + 100;
+    raw(1, 0, 1) = (1 << 14) + 200;
+    raw(2, 0, 1) = (1 << 14) + 300;
+
+    raw(0, 1, 0) = (1 << 14) + 37;
+    raw(1, 1, 0) = 38;
+    raw(2, 1, 0) = (3 << 14) + 39;
+
+    // gain 2
+    raw(0, 1, 1) = (3 << 14) + 100;
+    raw(1, 1, 1) = (3 << 14) + 200;
+    raw(2, 1, 1) = (3 << 14) + 300;
+
+    auto pedestal = calculate_pedestal<double>(raw.view(), 4);
+
+    REQUIRE(pedestal.size() == raw.size());
+    CHECK(pedestal(0, 0, 0) == 200);
+    CHECK(pedestal(1, 0, 0) == 0);
+    CHECK(pedestal(1, 0, 1) == 200);
+
+    auto pedestal_gain0 = calculate_pedestal_g0<double>(raw.view(), 4);
+
+    REQUIRE(pedestal_gain0.size() == 4);
+    CHECK(pedestal_gain0(0, 0) == 200);
+    CHECK(pedestal_gain0(1, 0) == 38);
+}