diff --git a/core/CMakeLists.txt b/core/CMakeLists.txt
index cbdfad8..cd4ef42 100644
--- a/core/CMakeLists.txt
+++ b/core/CMakeLists.txt
@@ -19,6 +19,8 @@ if(AARE_TESTS)
     set(TestSources
             ${CMAKE_CURRENT_SOURCE_DIR}/test/defs.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/test/ProducerConsumerQueue.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/test/NDArray.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/test/NDView.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/test/CircularFifo.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/test/wrappers.test.cpp
 
diff --git a/core/include/aare/NDView.hpp b/core/include/aare/NDView.hpp
index a08396d..b582613 100644
--- a/core/include/aare/NDView.hpp
+++ b/core/include/aare/NDView.hpp
@@ -6,6 +6,7 @@
 #include <numeric>
 #include <vector>
 
+template <ssize_t Ndim> using Shape = std::array<ssize_t, Ndim>;
 
 template <ssize_t Dim = 0, typename Strides> ssize_t element_offset(const Strides &) { return 0; }
 
diff --git a/core/test/NDArray.test.cpp b/core/test/NDArray.test.cpp
new file mode 100644
index 0000000..d9fa19a
--- /dev/null
+++ b/core/test/NDArray.test.cpp
@@ -0,0 +1,401 @@
+#include <catch2/catch_test_macros.hpp>
+#include "aare/NDArray.hpp"
+#include <array>
+
+
+// using reuss::DataSpan;
+// using reuss::Shape;
+
+TEST_CASE("Initial size is zero if no size is specified")
+{
+    NDArray<double> a;
+    REQUIRE(a.size() == 0);
+    REQUIRE(a.shape() == Shape<2>{0,0});
+}
+
+
+TEST_CASE("Construct from a DataSpan"){
+    std::vector<int> some_data(9,42);
+    NDView<int, 2> view(some_data.data(), Shape<2>{3,3});
+
+    NDArray<int,2> image(view);
+
+    REQUIRE(image.shape() == view.shape());
+    REQUIRE(image.size() == view.size());
+    REQUIRE(image.data() != view.data());
+
+    for (int i = 0; i<image.size(); ++i){
+        REQUIRE(image(i) == view(i));
+    }
+
+    //Changing the image doesn't change the view
+    image = 43;
+    for (int i = 0; i<image.size(); ++i){
+        REQUIRE(image(i) != view(i));
+    }
+
+
+}
+
+TEST_CASE("1D image"){
+    std::array<ssize_t, 1>shape{{20}};
+    NDArray<short,1>img(shape,3);
+    REQUIRE(img.size()==20);
+    REQUIRE(img(5)==3);
+
+}
+
+TEST_CASE("Accessing a const object"){
+    const NDArray<double, 3> img({3,4,5}, 0);
+    REQUIRE(img(1,1,1)==0);
+    REQUIRE(img.size() == 3*4*5);
+    REQUIRE(img.shape() == Shape<3>{3,4,5});
+    REQUIRE(img.shape(0) == 3);
+    REQUIRE(img.shape(1) == 4);
+    REQUIRE(img.shape(2) == 5);
+
+}
+
+TEST_CASE("Indexing of a 2D image")
+{
+    std::array<ssize_t, 2>shape{{ 3, 7 }};
+    NDArray<long> img( shape, 5 );
+    for (int i = 0; i != img.size(); ++i) {
+        REQUIRE(img(i) == 5);
+    }
+
+    for (int i = 0; i != img.size(); ++i) {
+        img(i) = i;
+    }
+    REQUIRE(img(0, 0) == 0);
+    REQUIRE(img(0, 1) == 1);
+    REQUIRE(img(1, 0) == 7);
+}
+
+TEST_CASE("Indexing of a 3D image")
+{
+    NDArray<float, 3> img{ {{ 3, 4, 2 }}, 5.0f };
+    for (int i = 0; i != img.size(); ++i) {
+        REQUIRE(img(i) == 5.0f);
+    }
+
+    //Double check general properties
+    REQUIRE(img.size() == 3 * 4 * 2);
+
+    for (int i = 0; i != img.size(); ++i) {
+        img(i) = float(i);
+    }
+    REQUIRE(img(0, 0, 0) == 0);
+    REQUIRE(img(0, 0, 1) == 1);
+    REQUIRE(img(0, 1, 1) == 3);
+    REQUIRE(img(1, 2, 0) == 12);
+    REQUIRE(img(2, 3, 1) == 23);
+}
+
+TEST_CASE("Divide double by int"){
+    NDArray<double,1> a{{5}, 5};
+    NDArray<int,1> b{{5},5};
+    a /=b;
+    for (auto it : a){
+        REQUIRE(it == 1.0);
+    }
+
+}
+
+TEST_CASE("Elementwise multiplication of 3D image")
+{
+    std::array<ssize_t, 3>shape{3, 4, 2};
+    NDArray<double, 3> a{ shape };
+    NDArray<double, 3> b{ shape };
+    for (int i = 0; i != a.size(); ++i) {
+        a(i) = i;
+        b(i) = i;
+    }
+    auto c = a * b;
+    REQUIRE(c(0, 0, 0) == 0 * 0);
+    REQUIRE(c(0, 0, 1) == 1 * 1);
+    REQUIRE(c(0, 1, 1) == 3 * 3);
+    REQUIRE(c(1, 2, 0) == 12 * 12);
+    REQUIRE(c(2, 3, 1) == 23 * 23);
+}
+
+TEST_CASE("Compare two images")
+{
+    NDArray<int> a;
+    NDArray<int> b;
+    CHECK(a == b);
+
+    a = NDArray<int>{ { 5, 10 }, 0 };
+    CHECK(a != b);
+
+    b = NDArray<int>{ { 5, 10 }, 0 };
+    CHECK(a == b);
+
+    b(3, 3) = 7;
+    CHECK(a != b);
+}
+
+TEST_CASE("Size and shape matches")
+{
+    ssize_t w = 15;
+    ssize_t h = 75;
+    std::array<ssize_t, 2> shape{ w, h };
+    NDArray<double> a{ shape };
+    REQUIRE(a.size() == w * h);
+    REQUIRE(a.shape() == shape);
+}
+
+TEST_CASE("Initial value matches for all elements")
+{
+    double v = 4.35;
+    NDArray<double> a{ { 5, 5 }, v };
+    for (int i = 0; i < a.size(); ++i) {
+        REQUIRE(a(i) == v);
+    }
+}
+
+TEST_CASE("Data layout of 3D image, fast index last"){
+    NDArray<int,3> a{{3,3,3}, 0};
+    REQUIRE(a.size()== 27);
+    int* ptr = a.data();
+
+    for (int i = 0; i<9; ++i){
+        *ptr++ = 10+i;
+        REQUIRE(a(0,0,i) == 10+i);
+        REQUIRE(a(i) == 10+i);
+
+    }
+}
+
+TEST_CASE("Bitwise and on data"){
+
+    NDArray<uint16_t, 1> a({3}, 0);
+    uint16_t mask = 0x3FF;
+    a(0) = 16684;
+    a(1) = 33068;
+    a(2) = 52608;
+
+    a &= mask;
+
+    REQUIRE(a(0) == 300);
+    REQUIRE(a(1) == 300);
+    REQUIRE(a(2) == 384);
+
+}
+
+// TEST_CASE("Benchmarks")
+// {
+//     NDArray<double> img;
+//     std::array<ssize_t, 2> shape{ 512, 1024 };
+//     BENCHMARK("Allocate 500k double image")
+//     {
+//         NDArray<double>im{ shape };
+//     }
+//     BENCHMARK("Allocate 500k double image with initial value")
+//     {
+//         NDArray<double>im{ shape, 3.14 };
+//     }
+
+//     NDArray<double> a{ shape, 1.2 };
+//     NDArray<double> b{ shape, 53. };
+//     auto c = a + b;
+//     c      = a * b;
+//     BENCHMARK("Multiply two images")
+//     {
+//         c = a * b;
+//     }
+//     BENCHMARK("Divide two images")
+//     {
+//         c = a / b;
+//     }
+//     BENCHMARK("Add two images")
+//     {
+//         c = a + b;
+//     }
+//     BENCHMARK("Subtract two images")
+//     {
+//         c = a - b;
+//     }
+// }
+
+TEST_CASE("Elementwise operatios on images")
+{
+    std::array<ssize_t, 2> shape{ 5, 5 };
+    double a_val = 3.0;
+    double b_val = 8.0;
+
+    SECTION("Add two images")
+    {
+        NDArray<double> A(shape, a_val);
+        NDArray<double> B(shape, b_val);
+
+        auto C = A + B;
+
+        //Value of C matches
+        for (int i = 0; i < C.size(); ++i) {
+            REQUIRE(C(i) == a_val + b_val);
+        }
+
+        //Value of A is not changed
+        for (int i = 0; i < A.size(); ++i) {
+            REQUIRE(A(i) == a_val);
+        }
+
+        //Value of B is not changed
+        for (int i = 0; i < B.size(); ++i) {
+            REQUIRE(B(i) == b_val);
+        }
+
+        //A, B and C referes to different data
+        REQUIRE(A.data() != B.data());
+        REQUIRE(B.data() != C.data());
+    }
+    SECTION("Subtract two images")
+    {
+        NDArray<double> A(shape, a_val);
+        NDArray<double> B(shape, b_val);
+        auto C = A - B;
+
+        //Value of C matches
+        for (int i = 0; i < C.size(); ++i) {
+            REQUIRE(C(i) == a_val - b_val);
+        }
+
+        //Value of A is not changed
+        for (int i = 0; i < A.size(); ++i) {
+            REQUIRE(A(i) == a_val);
+        }
+
+        //Value of B is not changed
+        for (int i = 0; i < B.size(); ++i) {
+            REQUIRE(B(i) == b_val);
+        }
+
+        //A, B and C referes to different data
+        REQUIRE(A.data() != B.data());
+        REQUIRE(B.data() != C.data());
+    }
+    SECTION("Multiply two images")
+    {
+        NDArray<double> A(shape, a_val);
+        NDArray<double> B(shape, b_val);
+        auto C = A * B;
+
+        //Value of C matches
+        for (int i = 0; i < C.size(); ++i) {
+            REQUIRE(C(i) == a_val * b_val);
+        }
+
+        //Value of A is not changed
+        for (int i = 0; i < A.size(); ++i) {
+            REQUIRE(A(i) == a_val);
+        }
+
+        //Value of B is not changed
+        for (int i = 0; i < B.size(); ++i) {
+            REQUIRE(B(i) == b_val);
+        }
+
+        //A, B and C referes to different data
+        REQUIRE(A.data() != B.data());
+        REQUIRE(B.data() != C.data());
+    }
+    SECTION("Divide two images")
+    {
+        NDArray<double> A(shape, a_val);
+        NDArray<double> B(shape, b_val);
+        auto C = A / B;
+
+        //Value of C matches
+        for (int i = 0; i < C.size(); ++i) {
+            REQUIRE(C(i) == a_val / b_val);
+        }
+
+        //Value of A is not changed
+        for (int i = 0; i < A.size(); ++i) {
+            REQUIRE(A(i) == a_val);
+        }
+
+        //Value of B is not changed
+        for (int i = 0; i < B.size(); ++i) {
+            REQUIRE(B(i) == b_val);
+        }
+
+        //A, B and C referes to different data
+        REQUIRE(A.data() != B.data());
+        REQUIRE(B.data() != C.data());
+    }
+
+    SECTION("subtract scalar")
+    {
+        NDArray<double> A(shape, a_val);
+        NDArray<double> B(shape, b_val);
+        double v = 1.0;
+        auto C   = A - v;
+        REQUIRE(C.data() != A.data());
+
+        //Value of C matches
+        for (int i = 0; i < C.size(); ++i) {
+            REQUIRE(C(i) == a_val - v);
+        }
+
+        //Value of A is not changed
+        for (int i = 0; i < A.size(); ++i) {
+            REQUIRE(A(i) == a_val);
+        }
+    }
+    SECTION("add scalar")
+    {
+        NDArray<double> A(shape, a_val);
+        NDArray<double> B(shape, b_val);
+        double v = 1.0;
+        auto C   = A + v;
+        REQUIRE(C.data() != A.data());
+
+        //Value of C matches
+        for (int i = 0; i < C.size(); ++i) {
+            REQUIRE(C(i) == a_val + v);
+        }
+
+        //Value of A is not changed
+        for (int i = 0; i < A.size(); ++i) {
+            REQUIRE(A(i) == a_val);
+        }
+    }
+    SECTION("divide with scalar")
+    {
+        NDArray<double> A(shape, a_val);
+        NDArray<double> B(shape, b_val);
+        double v = 3.7;
+        auto C   = A / v;
+        REQUIRE(C.data() != A.data());
+
+        //Value of C matches
+        for (int i = 0; i < C.size(); ++i) {
+            REQUIRE(C(i) == a_val / v);
+        }
+
+        //Value of A is not changed
+        for (int i = 0; i < A.size(); ++i) {
+            REQUIRE(A(i) == a_val);
+        }
+    }
+    SECTION("multiply with scalar")
+    {
+        NDArray<double> A(shape, a_val);
+        NDArray<double> B(shape, b_val);
+        double v = 3.7;
+        auto C   = A / v;
+        REQUIRE(C.data() != A.data());
+
+        //Value of C matches
+        for (int i = 0; i < C.size(); ++i) {
+            REQUIRE(C(i) == a_val / v);
+        }
+
+        //Value of A is not changed
+        for (int i = 0; i < A.size(); ++i) {
+            REQUIRE(A(i) == a_val);
+        }
+    }
+}
\ No newline at end of file
diff --git a/core/test/NDView.test.cpp b/core/test/NDView.test.cpp
new file mode 100644
index 0000000..9ea375f
--- /dev/null
+++ b/core/test/NDView.test.cpp
@@ -0,0 +1,178 @@
+#include "aare/NDView.hpp"
+#include <catch2/catch_test_macros.hpp>
+
+#include <iostream>
+#include <vector>
+
+
+
+TEST_CASE("Element reference 1D") {
+    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()));
+    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;
+    for (int i = 0; i != 12; ++i) {
+        vec.push_back(i);
+    }
+
+    NDView<int, 2> data(vec.data(), Shape<2>{3, 4});
+    REQUIRE(vec.size() == static_cast<size_t>(data.size()));
+    int i = 0;
+    for (int row = 0; row != 3; ++row) {
+        for (int col = 0; col != 4; ++col) {
+            REQUIRE(data(row, col) == i);
+            REQUIRE(data[i] == vec[i]);
+            ++i;
+        }
+    }
+}
+
+TEST_CASE("Element reference 3D") {
+    std::vector<int> vec;
+    for (int i = 0; i != 24; ++i) {
+        vec.push_back(i);
+    }
+    NDView<int, 3> data(vec.data(), Shape<3>{2, 3, 4});
+    REQUIRE(vec.size() == static_cast<size_t>(data.size()));
+    int i = 0;
+    for (int frame = 0; frame != 2; ++frame) {
+        for (int row = 0; row != 3; ++row) {
+            for (int col = 0; col != 4; ++col) {
+                REQUIRE(data(frame, row, col) == i);
+                REQUIRE(data[i] == vec[i]);
+                ++i;
+            }
+        }
+    }
+}
+
+TEST_CASE("Plus and miuns with single value") {
+    std::vector<int> vec;
+    for (int i = 0; i != 12; ++i) {
+        vec.push_back(i);
+    }
+    NDView<int, 2> data(vec.data(), Shape<2>{3, 4});
+    data += 5;
+    int i = 0;
+    for (int row = 0; row != 3; ++row) {
+        for (int col = 0; col != 4; ++col) {
+            REQUIRE(data(row, col) == i + 5);
+            ++i;
+        }
+    }
+    data -= 3;
+    i = 0;
+    for (int row = 0; row != 3; ++row) {
+        for (int col = 0; col != 4; ++col) {
+            REQUIRE(data(row, col) == i + 2);
+            ++i;
+        }
+    }
+}
+
+TEST_CASE("Multiply and divide with single value") {
+    std::vector<int> vec;
+    for (int i = 0; i != 12; ++i) {
+        vec.push_back(i);
+    }
+    NDView<int, 2> data(vec.data(), Shape<2>{3, 4});
+    data *= 5;
+    int i = 0;
+    for (int row = 0; row != 3; ++row) {
+        for (int col = 0; col != 4; ++col) {
+            REQUIRE(data(row, col) == i * 5);
+            ++i;
+        }
+    }
+    data /= 3;
+    i = 0;
+    for (int row = 0; row != 3; ++row) {
+        for (int col = 0; col != 4; ++col) {
+            REQUIRE(data(row, col) == (i * 5) / 3);
+            ++i;
+        }
+    }
+}
+
+TEST_CASE("elementwise assign"){
+    std::vector<int> vec(25);
+    NDView<int, 2> data(vec.data(), Shape<2>{5,5});
+
+    data = 3;
+    for (auto it : data){
+        REQUIRE(it == 3);
+    }
+}
+
+TEST_CASE("iterators") {
+    std::vector<int> vec;
+    for (int i = 0; i != 12; ++i) {
+        vec.push_back(i);
+    }
+    NDView<int, 1> data(vec.data(), Shape<1>{12});
+    int i = 0;
+    for (const auto item : data) {
+        REQUIRE(item == vec[i]);
+        ++i;
+    }
+    REQUIRE(i == 12);
+
+    for (auto ptr = data.begin(); ptr != data.end(); ++ptr) {
+        *ptr += 1;
+    }
+    for (auto& item : data){
+        ++item;
+    }
+
+    i = 0;
+    for (const auto item : data) {
+        REQUIRE(item == i + 2);
+        ++i;
+    }
+}
+
+
+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"){
+    std::vector<int> vec0{9,12,3};
+    std::vector<int> vec1{3,2,1};
+    std::vector<int> result{3,6,3};
+
+    NDView<int, 1> data0(vec0.data(), Shape<1>{static_cast<ssize_t>(vec0.size())});
+    NDView<int, 1> data1(vec1.data(), Shape<1>{static_cast<ssize_t>(vec1.size())});
+
+    data0 /= data1;
+    
+    for(size_t i =0; i!=vec0.size(); ++i){
+        REQUIRE(data0[i] == result[i] );
+    }
+}
+
+TEST_CASE("Retrieve shape"){
+    std::vector<int> vec;
+    for (int i = 0; i != 12; ++i) {
+        vec.push_back(i);
+    }
+    NDView<int,2> data(vec.data(), Shape<2>{3,4});
+    REQUIRE(data.shape()[0] == 3);
+    REQUIRE(data.shape()[1] == 4);
+
+}
\ No newline at end of file
diff --git a/core/test/ProducerConsumerQueue.test.cpp b/core/test/ProducerConsumerQueue.test.cpp
index a89ebad..03bba0c 100644
--- a/core/test/ProducerConsumerQueue.test.cpp
+++ b/core/test/ProducerConsumerQueue.test.cpp
@@ -42,4 +42,8 @@ TEST_CASE("Cannot pop from an empty queue"){
     folly::ProducerConsumerQueue<int> q(2);
     int a=0; 
     CHECK_FALSE(q.read(a));
-}
\ No newline at end of file
+}
+
+// TEST_CASE("fail"){
+//     CHECK(false);
+// }
\ No newline at end of file