updated release notes

- added ClusterVector::empty() to check if the vector is empty
- Fixed generation of missing docs for ClusterVector

CMakeLists.txt

@@ -369,22 +369,22 @@ set(PUBLICHEADERS
 
 set(SourceFiles
     ${CMAKE_CURRENT_SOURCE_DIR}/src/CtbRawFile.cpp
-    ${CMAKE_CURRENT_SOURCE_DIR}/src/defs.cpp
-    ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/decode.cpp
-    ${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/defs.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/DetectorGeometry.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/File.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/FilePtr.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/Fit.cpp
-    ${CMAKE_CURRENT_SOURCE_DIR}/src/DetectorGeometry.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/Interpolator.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/JungfrauDataFile.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.cpp
-    ${CMAKE_CURRENT_SOURCE_DIR}/src/Interpolator.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/PixelMap.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.cpp
-    ${CMAKE_CURRENT_SOURCE_DIR}/src/RawSubFile.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/RawSubFile.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/utils/task.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/utils/ifstream_helpers.cpp
 )

RELEASE.md

@@ -6,7 +6,11 @@
 Features:
 
 - Cluster finder now works with 5x5, 7x7 and 9x9 clusters
+- Added ClusterVector::empty() member
 
+Bugfixes:
+- Fixed reading RawFiles with ROI fully excluding some sub files. 
+- Removed unused file: ClusterFile.cpp 
 
 ### 2025.05.22
 
@@ -20,3 +24,4 @@ Bugfixes:
 
 
 
+

docs/CMakeLists.txt

@@ -45,12 +45,3 @@ add_custom_target(
     COMMENT "Generating documentation with Sphinx"
 )
 
-add_custom_target(
-    rst
-    COMMAND ${SPHINX_EXECUTABLE} -a -b html
-    -Dbreathe_projects.aare=${CMAKE_CURRENT_BINARY_DIR}/xml
-    -c "${SPHINX_BUILD}"
-    ${SPHINX_BUILD}/src 
-    ${SPHINX_BUILD}/html
-    COMMENT "Generating documentation with Sphinx"
-)

docs/Doxyfile.in

@@ -886,7 +886,7 @@ EXCLUDE_SYMLINKS       = NO
 # Note that the wildcards are matched against the file with absolute path, so to
 # exclude all test directories for example use the pattern */test/*
 
-EXCLUDE_PATTERNS       = */docs/* */tests/* */python/* */manual */slsDetectorServers/* */libs/* */integrationTests *README* */slsDetectorGui/* */ctbGui/* */slsDetectorCalibration/* *TobiSchluter*
+EXCLUDE_PATTERNS       = *build* */docs/* */tests/* *.test.cpp* */python/* */manual */slsDetectorServers/* */libs/* */integrationTests *README* *_deps* *TobiSchluter*
 
 # The EXCLUDE_SYMBOLS tag can be used to specify one or more symbol names
 # (namespaces, classes, functions, etc.) that should be excluded from the

docs/src/ClusterVector.rst

@@ -3,4 +3,13 @@ ClusterVector
 
 .. doxygenclass:: aare::ClusterVector
    :members:
-   :undoc-members:
+   :undoc-members:
+   :private-members:
+
+
+
+.. doxygenclass:: aare::ClusterVector< Cluster< T, ClusterSizeX, ClusterSizeY, CoordType > > 
+   :members:
+   :undoc-members:
+   :private-members:
+   

docs/src/Workflow.rst

@@ -37,7 +37,7 @@ unfamiliar steps.
 
 
 Checklists for deployment
-~~~~~~~~~~~~~~~~~~
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
 **Feature:**
 

include/aare/ClusterFile.hpp

@@ -14,7 +14,7 @@
 namespace aare {
 
 /*
-Binary cluster file. Expects data to be layed out as:
+Binary cluster file. Expects data to be laid 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

include/aare/ClusterVector.hpp

@@ -32,7 +32,8 @@ class ClusterVector; // Forward declaration
  */
 template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
           typename CoordType>
-class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
+class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> 
+{
 
     std::vector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> m_data{};
     int32_t m_frame_number{0}; // TODO! Check frame number size and type
@@ -122,6 +123,11 @@ class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
      */
     size_t size() const { return m_data.size(); }
 
+    /**
+     * @brief Check if the vector is empty
+     */
+    bool empty() const { return m_data.empty(); }
+
     uint8_t cluster_size_x() const { return ClusterSizeX; }
 
     uint8_t cluster_size_y() const { return ClusterSizeY; }

include/aare/NDArray.hpp

@@ -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.
@@ -380,7 +380,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();
@@ -423,7 +422,7 @@ template <typename T, ssize_t Ndim>
 void save(NDArray<T, Ndim> &img, std::string &pathname) {
     std::ofstream f;
     f.open(pathname, std::ios::binary);
-    f.write(reinterpret_cast<char *>(img.buffer()), img.size() * sizeof(T));
+    f.write(img.buffer(), img.size() * sizeof(T));
     f.close();
 }
 
@@ -433,39 +432,9 @@ NDArray<T, Ndim> load(const std::string &pathname,
     NDArray<T, Ndim> img{shape};
     std::ifstream f;
     f.open(pathname, std::ios::binary);
-    f.read(reinterpret_cast<char *>(img.buffer()), img.size() * sizeof(T));
+    f.read(img.buffer(), img.size() * sizeof(T));
     f.close();
     return img;
 }
 
-template <typename T, ssize_t Ndim = 1>
-NDArray<T, Ndim> load_non_binary_file(const std::string &filename,
-                                      const std::array<ssize_t, Ndim> shape) {
-    std::string word;
-    NDArray<T, Ndim> array(shape);
-
-    try {
-        std::ifstream file(filename, std::ios_base::in);
-        if (!file.good()) {
-            throw std::logic_error("file does not exist");
-        }
-
-        std::stringstream file_buffer;
-        file_buffer << file.rdbuf();
-
-        ssize_t counter = 0;
-        while (file_buffer >> word && counter < array.size()) {
-            array[counter] = static_cast<T>(
-                std::stod(word)); // TODO change for different Types
-            ++counter;
-        }
-
-        file.close();
-    } catch (const std::exception &e) {
-        std::cerr << "Error: " << e.what() << std::endl;
-    }
-
-    return array;
-}
-
 } // namespace aare

include/aare/NDView.hpp

@@ -6,7 +6,6 @@
 #include <array>
 #include <cassert>
 #include <cstdint>
-#include <fstream>
 #include <functional>
 #include <iomanip>
 #include <iostream>
@@ -68,20 +67,13 @@ class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
           size_(std::accumulate(std::begin(shape), std::end(shape), 1,
                                 std::multiplies<>())) {}
 
-    // NDView(T *buffer, const std::vector<ssize_t> &shape)
-    //     : buffer_(buffer),
-    //     strides_(c_strides<Ndim>(make_array<Ndim>(shape))),
-    //     shape_(make_array<Ndim>(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>
-    std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) const {
+    const std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) const {
         return buffer_[element_offset(strides_, index...)];
     }
 
@@ -93,33 +85,23 @@ 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) const { return buffer_[i]; }
-    T &operator[](ssize_t i) const { return buffer_[i]; }
+    T &operator()(ssize_t i)  { return buffer_[i]; }
+    T &operator[](ssize_t i)  { return buffer_[i]; }
+    const T &operator()(ssize_t i) const { return buffer_[i]; }
+    const T &operator[](ssize_t i) const { return buffer_[i]; }
 
     bool operator==(const NDView &other) const {
         if (size_ != other.size_)
             return false;
-        for (uint64_t i = 0; i != size_; ++i) {
+        if (shape_ != other.shape_)
+            return false;
+        for (size_t i = 0; i != size_; ++i) {
             if (buffer_[i] != other.buffer_[i])
                 return false;
         }
         return true;
     }
 
-    bool equals(const NDView<T, Ndim> &other, const T tolerance) const {
-        if (shape_ != other.shape_)
-            return false;
-
-        using SignedT = typename make_signed<T>::type;
-
-        for (uint32_t i = 0; i != size_; ++i)
-            if (std::abs(static_cast<SignedT>(buffer_[i]) -
-                         static_cast<SignedT>(other.buffer_[i])) > tolerance)
-                return false;
-
-        return true;
-    }
-
     NDView &operator+=(const T val) { return elemenwise(val, std::plus<T>()); }
     NDView &operator-=(const T val) { return elemenwise(val, std::minus<T>()); }
     NDView &operator*=(const T val) {
@@ -172,6 +154,7 @@ class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
     auto shape(ssize_t i) const { return shape_[i]; }
 
     T *data() { return buffer_; }
+    const T *data() const { return buffer_; }
     void print_all() const;
 
   private:
@@ -195,6 +178,7 @@ class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
         return *this;
     }
 };
+
 template <typename T, ssize_t Ndim> void NDView<T, Ndim>::print_all() const {
     for (auto row = 0; row < shape_[0]; ++row) {
         for (auto col = 0; col < shape_[1]; ++col) {
@@ -221,12 +205,4 @@ template <typename T> NDView<T, 1> make_view(std::vector<T> &vec) {
     return NDView<T, 1>(vec.data(), {static_cast<ssize_t>(vec.size())});
 }
 
-template <typename T, ssize_t Ndim>
-void save(NDView<T, Ndim> img, const std::string &pathname) {
-    std::ofstream f;
-    f.open(pathname, std::ios::binary);
-    f.write(reinterpret_cast<char *>(img.data()), img.size() * sizeof(T));
-    f.close();
-}
-
 } // namespace aare

include/aare/defs.hpp

@@ -9,7 +9,6 @@
 #include <stdexcept>
 #include <string>
 #include <string_view>
-#include <type_traits>
 #include <variant>
 #include <vector>
 
@@ -232,12 +231,4 @@ template <> FrameDiscardPolicy StringTo(const std::string & /*mode*/);
 
 using DataTypeVariants = std::variant<uint16_t, uint32_t>;
 
-template <typename T, bool = std::is_integral_v<T>> struct make_signed {
-    using type = T;
-};
-
-template <typename T> struct make_signed<T, true> {
-    using type = std::make_signed_t<T>;
-};
-
 } // namespace aare

src/ClusterFile.cpp

@@ -1,395 +0,0 @@
-#include "aare/ClusterFile.hpp"
-
-#include <algorithm>
-
-namespace aare {
-
-ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size,
-                         const std::string &mode)
-    : m_chunk_size(chunk_size), m_mode(mode) {
-
-    if (mode == "r") {
-        fp = fopen(fname.c_str(), "rb");
-        if (!fp) {
-            throw std::runtime_error("Could not open file for reading: " +
-                                     fname.string());
-        }
-    } else if (mode == "w") {
-        fp = fopen(fname.c_str(), "wb");
-        if (!fp) {
-            throw std::runtime_error("Could not open file for writing: " +
-                                     fname.string());
-        }
-    } else if (mode == "a") {
-        fp = fopen(fname.c_str(), "ab");
-        if (!fp) {
-            throw std::runtime_error("Could not open file for appending: " +
-                                     fname.string());
-        }
-    } else {
-        throw std::runtime_error("Unsupported mode: " + mode);
-    }
-}
-
-void ClusterFile::set_roi(ROI roi) { m_roi = roi; }
-
-void ClusterFile::set_noise_map(const NDView<int32_t, 2> noise_map) {
-    m_noise_map = NDArray<int32_t, 2>(noise_map);
-}
-
-void ClusterFile::set_gain_map(const NDView<double, 2> gain_map) {
-    m_gain_map = NDArray<double, 2>(gain_map);
-
-    // Gain map is passed as ADU/keV to avoid dividing in when applying the gain
-    // map we invert it here
-    for (auto &item : m_gain_map->view()) {
-        item = 1.0 / item;
-    }
-}
-
-ClusterFile::~ClusterFile() { close(); }
-
-void ClusterFile::close() {
-    if (fp) {
-        fclose(fp);
-        fp = nullptr;
-    }
-}
-
-void ClusterFile::write_frame(const ClusterVector<int32_t> &clusters) {
-    if (m_mode != "w" && m_mode != "a") {
-        throw std::runtime_error("File not opened for writing");
-    }
-    if (!(clusters.cluster_size_x() == 3) &&
-        !(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();
-    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();
-    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");
-    }
-}
-
-ClusterVector<int32_t> ClusterFile::read_clusters(size_t n_clusters) {
-    if (m_mode != "r") {
-        throw std::runtime_error("File not opened for reading");
-    }
-    if (m_noise_map || m_roi) {
-        return read_clusters_with_cut(n_clusters);
-    } else {
-        return read_clusters_without_cut(n_clusters);
-    }
-}
-
-ClusterVector<int32_t>
-ClusterFile::read_clusters_without_cut(size_t n_clusters) {
-    if (m_mode != "r") {
-        throw std::runtime_error("File not opened for reading");
-    }
-
-    ClusterVector<int32_t> clusters(3, 3, 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();
-    // 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;
-        }
-        nph_read += fread((buf + nph_read * clusters.item_size()),
-                          clusters.item_size(), nn, fp);
-        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)) {
-            clusters.set_frame_number(iframe);
-            // 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);
-                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);
-    if (m_gain_map)
-        clusters.apply_gain_map(m_gain_map->view());
-    return clusters;
-}
-
-ClusterVector<int32_t> ClusterFile::read_clusters_with_cut(size_t n_clusters) {
-    ClusterVector<int32_t> clusters(3, 3);
-    clusters.reserve(n_clusters);
-
-    // if there are photons left from previous frame read them first
-    if (m_num_left) {
-        while (m_num_left && clusters.size() < n_clusters) {
-            Cluster3x3 c = read_one_cluster();
-            if (is_selected(c)) {
-                clusters.push_back(c.x, c.y,
-                                   reinterpret_cast<std::byte *>(c.data));
-            }
-        }
-    }
-
-    // we did not have enough clusters left in the previous frame
-    // keep on reading frames until reaching n_clusters
-    if (clusters.size() < n_clusters) {
-        // sanity check
-        if (m_num_left) {
-            throw std::runtime_error(
-                LOCATION + "Entered second loop with clusters left\n");
-        }
-
-        int32_t frame_number = 0; // frame number needs to be 4 bytes!
-        while (fread(&frame_number, sizeof(frame_number), 1, fp)) {
-            if (fread(&m_num_left, sizeof(m_num_left), 1, fp)) {
-                clusters.set_frame_number(
-                    frame_number); // cluster vector will hold the last frame
-                                   // number
-                while (m_num_left && clusters.size() < n_clusters) {
-                    Cluster3x3 c = read_one_cluster();
-                    if (is_selected(c)) {
-                        clusters.push_back(
-                            c.x, c.y, reinterpret_cast<std::byte *>(c.data));
-                    }
-                }
-            }
-
-            // we have enough clusters, break out of the outer while loop
-            if (clusters.size() >= n_clusters)
-                break;
-        }
-    }
-    if (m_gain_map)
-        clusters.apply_gain_map(m_gain_map->view());
-
-    return clusters;
-}
-
-Cluster3x3 ClusterFile::read_one_cluster() {
-    Cluster3x3 c;
-    auto rc = fread(&c, sizeof(c), 1, fp);
-    if (rc != 1) {
-        throw std::runtime_error(LOCATION + "Could not read cluster");
-    }
-    --m_num_left;
-    return c;
-}
-
-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");
-    }
-    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("Could not read frame number");
-    }
-
-    if (fread(&m_num_left, sizeof(m_num_left), 1, fp) != 1) {
-        throw std::runtime_error("Could not read number of clusters");
-    }
-
-    ClusterVector<int32_t> clusters(3, 3);
-    clusters.reserve(m_num_left);
-    clusters.set_frame_number(frame_number);
-    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));
-        }
-    }
-    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) {
-        if (!(m_roi->contains(cl.x, cl.y))) {
-            return false;
-        }
-    }
-    if (m_noise_map) {
-        int32_t sum_1x1 = cl.data[4];   // central pixel
-        int32_t sum_2x2 = cl.sum_2x2(); // highest sum of 2x2 subclusters
-        int32_t sum_3x3 = cl.sum();     // sum of all pixels
-
-        auto noise =
-            (*m_noise_map)(cl.y, cl.x); // TODO! check if this is correct
-        if (sum_1x1 <= noise || sum_2x2 <= 2 * noise || sum_3x3 <= 3 * noise) {
-            return false;
-        }
-    }
-    // we passed all checks
-    return true;
-}
-
-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});
-
-    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;
-}
-
-/**
- * @brief Calculate the eta2 values for a 3x3 cluster and return them in a Eta2
- * struct containing etay, etax and the corner of the cluster.
- */
-Eta2 calculate_eta2(Cluster3x3 &cl) {
-    Eta2 eta{};
-
-    std::array<int32_t, 4> tot2;
-    tot2[0] = cl.data[0] + cl.data[1] + cl.data[3] + cl.data[4];
-    tot2[1] = cl.data[1] + cl.data[2] + cl.data[4] + cl.data[5];
-    tot2[2] = cl.data[3] + cl.data[4] + cl.data[6] + cl.data[7];
-    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)
-            eta.x = static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
-        if ((cl.data[1] + cl.data[4]) != 0)
-            eta.y = static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
-        eta.c = cBottomLeft;
-        break;
-    case cBottomRight:
-        if ((cl.data[2] + cl.data[5]) != 0)
-            eta.x = static_cast<double>(cl.data[5]) / (cl.data[4] + cl.data[5]);
-        if ((cl.data[1] + cl.data[4]) != 0)
-            eta.y = static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
-        eta.c = cBottomRight;
-        break;
-    case cTopLeft:
-        if ((cl.data[7] + cl.data[4]) != 0)
-            eta.x = static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
-        if ((cl.data[7] + cl.data[4]) != 0)
-            eta.y = static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
-        eta.c = cTopLeft;
-        break;
-    case cTopRight:
-        if ((cl.data[5] + cl.data[4]) != 0)
-            eta.x = static_cast<double>(cl.data[5]) / (cl.data[5] + cl.data[4]);
-        if ((cl.data[7] + cl.data[4]) != 0)
-            eta.y = static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
-        eta.c = cTopRight;
-        break;
-        // no default to allow compiler to warn about missing cases
-    }
-    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;
-}
-
-} // namespace aare

src/ClusterVector.test.cpp

@@ -7,9 +7,21 @@
 
 using aare::Cluster;
 using aare::ClusterVector;
+using C1 = Cluster<int32_t, 2, 2>;
+
+
+TEST_CASE("A newly created ClusterVector is empty") {
+    ClusterVector<C1> cv(4);
+    REQUIRE(cv.empty());
+}
+
+TEST_CASE("After pushing back one element the ClusterVector is not empty") {
+    ClusterVector<C1> cv(4);
+    cv.push_back(C1{1, 2, {3, 4}});
+    REQUIRE(!cv.empty());
+}
 
 TEST_CASE("item_size return the size of the cluster stored") {
-    using C1 = Cluster<int32_t, 2, 2>;
     ClusterVector<C1> cv(4);
     CHECK(cv.item_size() == sizeof(C1));
 
@@ -43,8 +55,7 @@ TEST_CASE("item_size return the size of the cluster stored") {
     CHECK(cv7.item_size() == sizeof(C7));
 }
 
-TEST_CASE("ClusterVector 2x2 int32_t capacity 4, push back then read",
-          "[.ClusterVector]") {
+TEST_CASE("ClusterVector 2x2 int32_t capacity 4, push back then read") {
 
     ClusterVector<Cluster<int32_t, 2, 2>> cv(4);
     REQUIRE(cv.capacity() == 4);
@@ -70,7 +81,7 @@ TEST_CASE("ClusterVector 2x2 int32_t capacity 4, push back then read",
     }
 }
 
-TEST_CASE("Summing 3x1 clusters of int64", "[.ClusterVector]") {
+TEST_CASE("Summing 3x1 clusters of int64") {
     ClusterVector<Cluster<int32_t, 3, 1>> cv(2);
     REQUIRE(cv.capacity() == 2);
     REQUIRE(cv.size() == 0);
@@ -102,7 +113,7 @@ TEST_CASE("Summing 3x1 clusters of int64", "[.ClusterVector]") {
     */
 }
 
-TEST_CASE("Storing floats", "[.ClusterVector]") {
+TEST_CASE("Storing floats") {
     ClusterVector<Cluster<float, 2, 4>> cv(10);
     REQUIRE(cv.capacity() == 10);
     REQUIRE(cv.size() == 0);
@@ -129,7 +140,7 @@ TEST_CASE("Storing floats", "[.ClusterVector]") {
     */
 }
 
-TEST_CASE("Push back more than initial capacity", "[.ClusterVector]") {
+TEST_CASE("Push back more than initial capacity") {
 
     ClusterVector<Cluster<int32_t, 2, 2>> cv(2);
     auto initial_data = cv.data();
@@ -162,8 +173,7 @@ TEST_CASE("Push back more than initial capacity", "[.ClusterVector]") {
     REQUIRE(initial_data != cv.data());
 }
 
-TEST_CASE("Concatenate two cluster vectors where the first has enough capacity",
-          "[.ClusterVector]") {
+TEST_CASE("Concatenate two cluster vectors where the first has enough capacity") {
     ClusterVector<Cluster<int32_t, 2, 2>> cv1(12);
     Cluster<int32_t, 2, 2> c1 = {1, 2, {3, 4, 5, 6}};
     cv1.push_back(c1);
@@ -192,8 +202,7 @@ TEST_CASE("Concatenate two cluster vectors where the first has enough capacity",
     REQUIRE(ptr[3].y == 17);
 }
 
-TEST_CASE("Concatenate two cluster vectors where we need to allocate",
-          "[.ClusterVector]") {
+TEST_CASE("Concatenate two cluster vectors where we need to allocate") {
     ClusterVector<Cluster<int32_t, 2, 2>> cv1(2);
     Cluster<int32_t, 2, 2> c1 = {1, 2, {3, 4, 5, 6}};
     cv1.push_back(c1);
@@ -229,7 +238,7 @@ struct ClusterTestData {
     std::vector<int64_t> index_map_y;
 };
 
-TEST_CASE("Gain Map Calculation Index Map", "[.ClusterVector][.gain_map]") {
+TEST_CASE("Gain Map Calculation Index Map") {
 
     auto clustertestdata = GENERATE(
         ClusterTestData{3,

src/NDView.test.cpp

@@ -21,6 +21,57 @@ TEST_CASE("Element reference 1D") {
     }
 }
 
+
+TEST_CASE("Assign elements through () and []") {
+    std::vector<int> vec;
+    for (int i = 0; i != 10; ++i) {
+        vec.push_back(i);
+    }
+    NDView<int, 1> data(vec.data(), Shape<1>{10});
+    REQUIRE(vec.size() == static_cast<size_t>(data.size()));
+    
+    data[3] = 187;
+    data(4) = 512;
+
+
+    REQUIRE(data(0) == 0);
+    REQUIRE(data[0] == 0);
+    REQUIRE(data(1) == 1);
+    REQUIRE(data[1] == 1);
+    REQUIRE(data(2) == 2);
+    REQUIRE(data[2] == 2);
+    REQUIRE(data(3) == 187);
+    REQUIRE(data[3] == 187);
+    REQUIRE(data(4) == 512);
+    REQUIRE(data[4] == 512);
+    REQUIRE(data(5) == 5);
+    REQUIRE(data[5] == 5);
+    REQUIRE(data(6) == 6);
+    REQUIRE(data[6] == 6);
+    REQUIRE(data(7) == 7);
+    REQUIRE(data[7] == 7);
+    REQUIRE(data(8) == 8);
+    REQUIRE(data[8] == 8);
+    REQUIRE(data(9) == 9);
+    REQUIRE(data[9] == 9);
+
+
+}
+
+TEST_CASE("Element reference 1D with a const NDView") {
+    std::vector<int> vec;
+    for (int i = 0; i != 10; ++i) {
+        vec.push_back(i);
+    }
+    const NDView<int, 1> data(vec.data(), Shape<1>{10});
+    REQUIRE(vec.size() == static_cast<size_t>(data.size()));
+    for (int i = 0; i != 10; ++i) {
+        REQUIRE(data(i) == vec[i]);
+        REQUIRE(data[i] == vec[i]);
+    }
+}
+
+
 TEST_CASE("Element reference 2D") {
     std::vector<int> vec(12);
     std::iota(vec.begin(), vec.end(), 0);
@@ -56,7 +107,7 @@ TEST_CASE("Element reference 3D") {
     }
 }
 
-TEST_CASE("Plus and miuns with single value") {
+TEST_CASE("Plus and minus with single value") {
     std::vector<int> vec(12);
     std::iota(vec.begin(), vec.end(), 0);
     NDView<int, 2> data(vec.data(), Shape<2>{3, 4});
@@ -137,16 +188,9 @@ TEST_CASE("iterators") {
     }
 }
 
-// TEST_CASE("shape from vector") {
-//     std::vector<int> vec;
-//     for (int i = 0; i != 12; ++i) {
-//         vec.push_back(i);
-//     }
-//     std::vector<ssize_t> shape{3, 4};
-//     NDView<int, 2> data(vec.data(), shape);
-// }
 
-TEST_CASE("divide with another span") {
+
+TEST_CASE("divide with another NDView") {
     std::vector<int> vec0{9, 12, 3};
     std::vector<int> vec1{3, 2, 1};
     std::vector<int> result{3, 6, 3};
@@ -183,6 +227,30 @@ TEST_CASE("compare two views") {
     REQUIRE((view1 == view2));
 }
 
+TEST_CASE("Compare two views with different size"){
+    std::vector<int> vec1(12);
+    std::iota(vec1.begin(), vec1.end(), 0);
+    NDView<int, 2> view1(vec1.data(), Shape<2>{3, 4});
+
+    std::vector<int> vec2(8);
+    std::iota(vec2.begin(), vec2.end(), 0);
+    NDView<int, 2> view2(vec2.data(), Shape<2>{2, 4});
+
+    REQUIRE_FALSE(view1 == view2);
+}
+
+TEST_CASE("Compare two views with same size but different shape"){
+    std::vector<int> vec1(12);
+    std::iota(vec1.begin(), vec1.end(), 0);
+    NDView<int, 2> view1(vec1.data(), Shape<2>{3, 4});
+
+    std::vector<int> vec2(12);
+    std::iota(vec2.begin(), vec2.end(), 0);
+    NDView<int, 2> view2(vec2.data(), Shape<2>{2, 6});
+
+    REQUIRE_FALSE(view1 == view2);
+}
+
 TEST_CASE("Create a view over a vector") {
     std::vector<int> vec(12);
     std::iota(vec.begin(), vec.end(), 0);