Merge branch 'main' into dev/reduce

include/aare/ClusterFinder.hpp

@@ -144,9 +144,9 @@ class ClusterFinder {
                                     static_cast<CT>(
                                         m_pedestal.mean(iy + ir, ix + ic));
                                 cluster.data[i] =
-                                    tmp; // Watch for out of bounds access
-                                i++;
+                                    tmp; // Watch for out of bounds access                                
                             }
+                            i++;
                         }
                     }
 

include/aare/Frame.hpp

@@ -105,7 +105,7 @@ class Frame {
      * @tparam T type of the pixels
      * @return NDView<T, 2>
      */
-    template <typename T> NDView<T, 2> view() {
+    template <typename T> NDView<T, 2> view() & {
         std::array<ssize_t, 2> shape = {static_cast<ssize_t>(m_rows),
                                         static_cast<ssize_t>(m_cols)};
         T *data = reinterpret_cast<T *>(m_data);

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,
-                                std::multiplies<>())),
+          size_(num_elements(shape_)),
           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_)),
@@ -380,12 +397,6 @@ NDArray<T, Ndim> NDArray<T, Ndim>::operator*(const T &value) {
     result *= value;
     return result;
 }
-// template <typename T, ssize_t Ndim> void NDArray<T, Ndim>::Print() {
-//     if (shape_[0] < 20 && shape_[1] < 20)
-//         Print_all();
-//     else
-//         Print_some();
-// }
 
 template <typename T, ssize_t Ndim>
 std::ostream &operator<<(std::ostream &os, const NDArray<T, Ndim> &arr) {
@@ -437,4 +448,23 @@ 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) {
+    if (numerator.shape() != denominator.shape()) {
+        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]);
+        } else {
+            result[i] = RT{0}; // or handle division by zero as needed
+        }
+    }
+    return result;
+}
+
 } // namespace aare

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;
@@ -66,17 +93,28 @@ class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
         : buffer_(buffer), strides_(c_strides<Ndim>(shape)), shape_(shape),
           size_(std::accumulate(std::begin(shape), std::end(shape), 1,
                                 std::multiplies<>())) {}
-
+                   
     template <typename... Ix>
     std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) {
         return buffer_[element_offset(strides_, index...)];
     }
 
     template <typename... Ix>
-    const std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) const {
+    std::enable_if_t<sizeof...(Ix) == 1 && (Ndim > 1), NDView<T, Ndim - 1>> operator()(Ix... index) {
+        // return a view of the next dimension
+        std::array<ssize_t, Ndim - 1> new_shape{};
+        std::copy_n(shape_.begin() + 1, Ndim - 1, new_shape.begin());
+        return NDView<T, Ndim - 1>(&buffer_[element_offset(strides_, index...)],
+                                   new_shape);
+        
+    }
+
+    template <typename... Ix>
+    std::enable_if_t<sizeof...(Ix) == Ndim, const T &> operator()(Ix... index) const {
         return buffer_[element_offset(strides_, index...)];
     }
 
+
     ssize_t size() const { return static_cast<ssize_t>(size_); }
     size_t total_bytes() const { return size_ * sizeof(T); }
     std::array<ssize_t, Ndim> strides() const noexcept { return strides_; }
@@ -85,9 +123,19 @@ class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
     T *end() { return buffer_ + size_; }
     T const *begin() const { return buffer_; }
     T const *end() const { return buffer_ + size_; }
-    T &operator()(ssize_t i)  { return buffer_[i]; }
+    
+    
+
+
+
+    /**
+     * @brief Access element at index i.
+     */
     T &operator[](ssize_t i)  { return buffer_[i]; }
-    const T &operator()(ssize_t i) const { return buffer_[i]; }
+
+    /**
+     * @brief Access element at index i.
+     */
     const T &operator[](ssize_t i) const { return buffer_[i]; }
 
     bool operator==(const NDView &other) const {
@@ -157,6 +205,22 @@ class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
     const T *data() const { return buffer_; }
     void print_all() const;
 
+    /**
+     * @brief Create a subview of a range of the first dimension. 
+     * This is useful for splitting a batches of frames in parallel processing.
+     * @param first The first index of the subview (inclusive).
+     * @param last The last index of the subview (exclusive).
+     * @return A new NDView that is a subview of the current view.
+     * @throws std::runtime_error if the range is invalid.
+     */
+    NDView sub_view(ssize_t first, ssize_t last) const {
+        if (first < 0 || last > shape_[0] || first >= last)
+            throw std::runtime_error(LOCATION + "Invalid sub_view range");
+        auto new_shape = shape_;
+        new_shape[0] = last - first;
+        return NDView(buffer_ + first * strides_[0], new_shape);
+    }
+
   private:
     T *buffer_{nullptr};
     std::array<ssize_t, Ndim> strides_{};

include/aare/VarClusterFinder.hpp

@@ -240,14 +240,14 @@ template <typename T> void VarClusterFinder<T>::first_pass() {
 
     for (ssize_t i = 0; i < original_.size(); ++i) {
         if (use_noise_map)
-            threshold_ = 5 * noiseMap(i);
-        binary_(i) = (original_(i) > threshold_);
+            threshold_ = 5 * noiseMap[i];
+        binary_[i] = (original_[i] > threshold_);
     }
 
     for (int i = 0; i < shape_[0]; ++i) {
         for (int j = 0; j < shape_[1]; ++j) {
 
-            // do we have someting to process?
+            // do we have something to process?
             if (binary_(i, j)) {
                 auto tmp = check_neighbours(i, j);
                 if (tmp != 0) {

include/aare/calibration.hpp

@@ -1,6 +1,9 @@
 #pragma once
 
+#include "aare/NDArray.hpp"
+#include "aare/NDView.hpp"
 #include "aare/defs.hpp"
+#include "aare/utils/par.hpp"
 #include "aare/utils/task.hpp"
 #include <cstdint>
 #include <future>
@@ -55,32 +58,152 @@ ALWAYS_INLINE std::pair<uint16_t, int16_t> get_value_and_gain(uint16_t raw) {
 
 template <class T>
 void apply_calibration_impl(NDView<T, 3> res, NDView<uint16_t, 3> raw_data,
-                       NDView<T, 3> ped, NDView<T, 3> cal, int start,
-                       int stop) {
+                            NDView<T, 3> ped, NDView<T, 3> cal, int start,
+                            int stop) {
 
     for (int frame_nr = start; frame_nr != stop; ++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));
+
+                // Using multiplication does not seem to speed up the code here
+                // ADU/keV is the standard unit for the calibration which
+                // means rewriting the formula is not worth it.
                 res(frame_nr, row, col) =
-                    (value - ped(gain, row, col)) / cal(gain, row, col); //TODO! use multiplication
+                    (value - ped(gain, row, col)) / cal(gain, row, col);
             }
         }
     }
 }
 
 template <class T>
+void apply_calibration_impl(NDView<T, 3> res, NDView<uint16_t, 3> raw_data,
+                            NDView<T, 2> ped, NDView<T, 2> cal, int start,
+                            int stop) {
+
+    for (int frame_nr = start; frame_nr != stop; ++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));
+
+                // Using multiplication does not seem to speed up the code here
+                // ADU/keV is the standard unit for the calibration which
+                // means rewriting the formula is not worth it.
+
+                // Set the value to 0 if the gain is not 0
+                if (gain == 0)
+                    res(frame_nr, row, col) =
+                        (value - ped(row, col)) / cal(row, col);
+                else
+                    res(frame_nr, row, col) = 0;
+            }
+        }
+    }
+}
+
+template <class T, ssize_t Ndim = 3>
 void apply_calibration(NDView<T, 3> res, NDView<uint16_t, 3> raw_data,
-                       NDView<T, 3> ped, NDView<T, 3> cal,
+                       NDView<T, Ndim> ped, NDView<T, Ndim> cal,
                        ssize_t n_threads = 4) {
     std::vector<std::future<void>> futures;
     futures.reserve(n_threads);
     auto limits = split_task(0, raw_data.shape(0), n_threads);
     for (const auto &lim : limits)
-        futures.push_back(std::async(&apply_calibration_impl<T>, res, raw_data, ped, cal,
-                                     lim.first, lim.second));
+        futures.push_back(std::async(
+            static_cast<void (*)(NDView<T, 3>, NDView<uint16_t, 3>,
+                                 NDView<T, Ndim>, NDView<T, Ndim>, 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) {
+    constexpr ssize_t num_gains = only_gain0 ? 1 : 3;
+    NDArray<size_t, 3> accumulator(
+        std::array<ssize_t, 3>{num_gains, raw_data.shape(1), raw_data.shape(2)},
+        0);
+    NDArray<size_t, 3> count(
+        std::array<ssize_t, 3>{num_gains, 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 && only_gain0)
+                    continue;
+                accumulator(gain, row, col) += value;
+                count(gain, row, col) += 1;
+            }
+        }
+    }
+
+    return {std::move(accumulator), std::move(count)};
+}
+
+template <typename T, bool only_gain0 = false>
+NDArray<T, 3 - static_cast<ssize_t>(only_gain0)>
+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);
+
+    auto subviews = make_subviews(raw_data, n_threads);
+
+    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<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);
+
+    // Combine the results from the futures
+    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);
+
+}
+
+/**
+ * @brief Count the number of switching pixels in the raw data.
+ * This function counts the number of pixels that switch between G1 and G2 gain.
+ * It returns an NDArray with the number of switching pixels per pixel.
+ * @param raw_data The NDView containing the raw data
+ * @return An NDArray with the number of switching pixels per pixel
+ */
+NDArray<int, 2> count_switching_pixels(NDView<uint16_t, 3> raw_data);
+
+/**
+ * @brief Count the number of switching pixels in the raw data.
+ * This function counts the number of pixels that switch between G1 and G2 gain.
+ * It returns an NDArray with the number of switching pixels per pixel.
+ * @param raw_data The NDView containing the raw data
+ * @param n_threads The number of threads to use for parallel processing
+ * @return An NDArray with the number of switching pixels per pixel
+ */
+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

include/aare/utils/par.hpp

@@ -1,7 +1,10 @@
+#pragma once
 #include <thread>
 #include <utility>
 #include <vector>
 
+#include "aare/utils/task.hpp"
+
 namespace aare {
 
 template <typename F>
@@ -15,4 +18,17 @@ void RunInParallel(F func, const std::vector<std::pair<int, int>> &tasks) {
         thread.join();
     }
 }
+
+
+template <typename T>
+std::vector<NDView<T,3>> make_subviews(NDView<T, 3> &data, ssize_t n_threads) {
+    std::vector<NDView<T, 3>> subviews;
+    subviews.reserve(n_threads);
+    auto limits = split_task(0, data.shape(0), n_threads);
+    for (const auto &lim : limits) {
+        subviews.push_back(data.sub_view(lim.first, lim.second));
+    }
+    return subviews;
+}
+
 } // namespace aare

include/aare/utils/task.hpp

@@ -1,4 +1,4 @@
-
+#pragma once
 #include <utility>
 #include <vector>