WIP

CMakeLists.txt

@@ -0,0 +1,367 @@
+cmake_minimum_required(VERSION 3.14)
+
+project(aare 
+    VERSION 1.0.0
+    DESCRIPTION "Data processing library for PSI detectors"
+    HOMEPAGE_URL "https://github.com/slsdetectorgroup/aare"
+    LANGUAGES C CXX 
+)
+
+set(CMAKE_CXX_STANDARD 17) 
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_CXX_EXTENSIONS OFF)
+
+if (${CMAKE_VERSION} VERSION_GREATER "3.24")
+    cmake_policy(SET CMP0135 NEW) #Fetch content download timestamp
+endif()
+cmake_policy(SET CMP0079 NEW)
+
+include(GNUInstallDirs)
+include(FetchContent)
+
+#Set default build type if none was specified
+include(cmake/helpers.cmake)
+default_build_type("Release")
+
+
+
+set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH})
+
+option(AARE_USE_WARNINGS "Enable warnings" ON)
+option(AARE_PYTHON_BINDINGS "Build python bindings" ON)
+option(AARE_TESTS "Build tests" ON)
+option(AARE_EXAMPLES "Build examples" ON)
+option(AARE_IN_GITHUB_ACTIONS "Running in Github Actions" OFF)
+option(AARE_DOCS "Build documentation" OFF)
+
+option(AARE_FETCH_FMT "Use FetchContent to download fmt" ON)
+option(AARE_FETCH_PYBIND11 "Use FetchContent to download pybind11" ON)
+option(AARE_FETCH_CATCH "Use FetchContent to download catch2" ON)
+option(AARE_FETCH_JSON "Use FetchContent to download nlohmann::json" ON)
+option(AARE_FETCH_ZMQ "Use FetchContent to download libzmq" ON)
+option(ENABLE_DRAFTS "Enable zmq drafts (depends on gnutls or nss)" OFF)
+
+
+#Convenience option to use system libraries
+option(AARE_SYSTEM_LIBRARIES "Use system libraries" OFF)
+if(AARE_SYSTEM_LIBRARIES)
+    message(STATUS "Build using system libraries")
+    set(AARE_FETCH_FMT OFF CACHE BOOL "Disabled FetchContent for FMT" FORCE)
+    set(AARE_FETCH_PYBIND11 OFF CACHE BOOL "Disabled FetchContent for pybind11" FORCE)
+    set(AARE_FETCH_CATCH OFF CACHE BOOL "Disabled FetchContent for catch2" FORCE)
+    set(AARE_FETCH_JSON OFF CACHE BOOL "Disabled FetchContent for nlohmann::json" FORCE)
+    set(AARE_FETCH_ZMQ OFF CACHE BOOL "Disabled FetchContent for libzmq" FORCE)
+
+endif()
+
+
+set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
+
+if(AARE_FETCH_ZMQ)
+    FetchContent_Declare(
+        libzmq
+        GIT_REPOSITORY https://github.com/zeromq/libzmq.git
+        GIT_TAG        v4.3.4
+    )
+    # Disable unwanted options from libzmq
+    set(BUILD_TESTS OFF CACHE BOOL "Switch off libzmq test build")
+    set(BUILD_SHARED OFF CACHE BOOL "Switch off libzmq shared libs")
+    set(WITH_PERF_TOOL OFF CACHE BOOL "")
+    set(ENABLE_CPACK OFF CACHE BOOL "")
+    set(ENABLE_CLANG OFF CACHE BOOL "")
+    set(ENABLE_CURVE OFF CACHE BOOL "")
+    set(ENABLE_DRAFTS OFF CACHE BOOL "")
+
+    # TODO! Verify that this is what we want to do in aare
+    # Using GetProperties and Populate to be able to exclude zmq
+    # from install (not possible with FetchContent_MakeAvailable(libzmq))
+    FetchContent_GetProperties(libzmq)
+    if(NOT libzmq_POPULATED)
+        FetchContent_Populate(libzmq)
+        add_subdirectory(${libzmq_SOURCE_DIR} ${libzmq_BINARY_DIR} EXCLUDE_FROM_ALL)
+    endif()
+else()
+    find_package(ZeroMQ 4 REQUIRED)
+endif()
+
+
+if (AARE_FETCH_FMT)
+    set(FMT_TEST OFF CACHE INTERNAL "disabling fmt tests")
+    FetchContent_Declare(
+            fmt
+            GIT_REPOSITORY  https://github.com/fmtlib/fmt.git
+            GIT_TAG         10.2.1
+            GIT_PROGRESS    TRUE
+            USES_TERMINAL_DOWNLOAD TRUE
+    )
+    FetchContent_MakeAvailable(fmt)
+    set_property(TARGET fmt PROPERTY POSITION_INDEPENDENT_CODE ON)
+else()
+      find_package(fmt 6 REQUIRED)
+endif()
+
+#TODO! Add options for nlohmann_json as well
+find_package(nlohmann_json 3.11.3 REQUIRED)
+
+
+include(GNUInstallDirs)
+
+# If conda build, always set lib dir to 'lib'
+if($ENV{CONDA_BUILD})
+  set(CMAKE_INSTALL_LIBDIR "lib")
+endif()
+
+# Set lower / upper case project names
+string(TOUPPER "${PROJECT_NAME}" PROJECT_NAME_UPPER)
+string(TOLOWER "${PROJECT_NAME}" PROJECT_NAME_LOWER)
+
+
+# Set targets export name (used by slsDetectorPackage and dependencies)
+set(TARGETS_EXPORT_NAME "${PROJECT_NAME_LOWER}-targets")
+set(namespace "aare::")
+
+set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH})
+
+# Check if project is being used directly or via add_subdirectory
+set(AARE_MASTER_PROJECT OFF)
+if (CMAKE_CURRENT_SOURCE_DIR STREQUAL CMAKE_SOURCE_DIR)
+    set(AARE_MASTER_PROJECT ON)
+endif()
+
+add_library(aare_compiler_flags INTERFACE)
+target_compile_features(aare_compiler_flags INTERFACE cxx_std_17)
+
+if(AARE_PYTHON_BINDINGS)
+    add_subdirectory(python)
+endif()
+
+#################
+# MSVC specific #
+#################
+if(MSVC)
+add_compile_definitions(AARE_MSVC)
+if(CMAKE_BUILD_TYPE STREQUAL "Release")
+    message(STATUS "Release build")
+    target_compile_options(aare_compiler_flags INTERFACE /O2)
+else()
+    message(STATUS "Debug build")
+    target_compile_options(
+        aare_compiler_flags 
+        INTERFACE 
+            /Od 
+            /Zi 
+            /MDd
+            /D_ITERATOR_DEBUG_LEVEL=2
+    )
+    target_link_options(
+        aare_compiler_flags
+        INTERFACE
+            /DEBUG:FULL
+            )
+endif()
+target_compile_options(
+    aare_compiler_flags 
+    INTERFACE
+    /w # disable warnings
+)
+
+
+else()
+######################
+# GCC/Clang specific #
+######################
+
+if(CMAKE_BUILD_TYPE STREQUAL "Release")
+    message(STATUS "Release build")
+    target_compile_options(aare_compiler_flags INTERFACE -O3)
+else()
+    message(STATUS "Debug build")
+    target_compile_options(
+        aare_compiler_flags 
+        INTERFACE 
+            -Og 
+            -ggdb3 
+            # -D_GLIBCXX_DEBUG # causes errors with boost
+            -D_GLIBCXX_DEBUG_PEDANTIC
+    )
+
+endif()
+
+if(AARE_USE_WARNINGS)
+        target_compile_options(
+            aare_compiler_flags 
+            INTERFACE 
+                -Wall 
+                -Wextra 
+                -pedantic 
+                -Wshadow 
+                -Wnon-virtual-dtor
+                -Woverloaded-virtual
+                -Wdouble-promotion
+                -Wformat=2
+                -Wredundant-decls
+                -Wvla
+                -Wdouble-promotion
+                -Werror=return-type #important can cause segfault in optimzed builds
+            )
+endif()
+
+
+
+endif() #GCC/Clang specific
+
+
+
+
+
+
+
+
+if(AARE_TESTS)
+    enable_testing()
+    add_subdirectory(tests)
+endif()
+
+###------------------------------------------------------------------------------MAIN LIBRARY
+###------------------------------------------------------------------------------------------
+
+set(PUBLICHEADERS
+    include/aare/ClusterFinder.hpp
+    include/aare/defs.hpp
+    include/aare/Dtype.hpp
+    include/aare/File.hpp
+    include/aare/FileInterface.hpp
+    include/aare/Frame.hpp
+    include/aare/json.hpp
+    include/aare/NDArray.hpp
+    include/aare/NDView.hpp
+    include/aare/NumpyFile.hpp
+    include/aare/NumpyHelpers.hpp
+    include/aare/Pedestal.hpp
+    include/aare/RawFile.hpp
+    include/aare/SubFile.hpp
+    include/aare/VarClusterFinder.hpp
+
+)
+
+
+set(SourceFiles
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/defs.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/File.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/SubFile.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.cpp
+)
+
+
+add_library(aare_core STATIC ${SourceFiles})
+target_include_directories(aare_core PUBLIC 
+    "$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>" 
+    "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
+)
+
+target_link_libraries(aare_core PUBLIC fmt::fmt PRIVATE aare_compiler_flags nlohmann_json::nlohmann_json)
+
+set_target_properties(aare_core PROPERTIES
+    # ARCHIVE_OUTPUT_NAME SlsDetectorStatic
+    ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}
+    PUBLIC_HEADER "${PUBLICHEADERS}"
+)
+
+if (AARE_PYTHON_BINDINGS)
+    set_property(TARGET aare_core PROPERTY POSITION_INDEPENDENT_CODE ON)
+endif()
+
+if(AARE_TESTS)
+    set(TestSources
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/defs.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/NDArray.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/NDView.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinder.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/Pedestal.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.test.cpp
+
+        )
+    target_sources(tests PRIVATE ${TestSources} )
+endif()
+
+
+
+###------------------------------------------------------------------------------------------
+###------------------------------------------------------------------------------------------
+
+
+if(AARE_MASTER_PROJECT)
+    install(TARGETS aare_core aare_compiler_flags 
+        EXPORT "${TARGETS_EXPORT_NAME}"
+        LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
+        ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
+        PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/aare
+    )
+endif()
+
+set(CMAKE_POSITION_INDEPENDENT_CODE ON)
+set(CMAKE_INSTALL_RPATH $ORIGIN)
+set(CMAKE_BUILD_WITH_INSTALL_RPATH FALSE)
+
+
+#Overall target to link to when using the library
+add_library(aare INTERFACE) 
+target_link_libraries(aare INTERFACE aare_core aare_compiler_flags)
+target_include_directories(aare INTERFACE 
+    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
+    $<INSTALL_INTERFACE:include>
+)
+
+# add_subdirectory(examples)
+
+if(AARE_DOCS)
+    add_subdirectory(docs)
+endif()
+
+
+# custom target to run check formatting with clang-format
+add_custom_target(
+    check-format
+    COMMAND find \( -name "*.cpp" -o -name "*.hpp" \) -not -path "./build/*" | xargs -I {} -n 1 -P 10 bash -c  "clang-format -Werror -style=\"file:.clang-format\" {} | diff {} -"
+    WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
+    COMMENT "Checking code formatting with clang-format"
+    VERBATIM
+    
+)
+
+add_custom_target(
+    format-files
+    COMMAND find \( -name "*.cpp" -o -name "*.hpp" \) -not -path "./build/*" | xargs -I {} -n 1 -P 10 bash -c "clang-format -i -style=\"file:.clang-format\" {}"
+    WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
+    COMMENT "Formatting with clang-format"
+    VERBATIM
+)
+
+if (AARE_IN_GITHUB_ACTIONS)
+    message(STATUS "Running in Github Actions")
+    set(CLANG_TIDY_COMMAND "clang-tidy-17")
+else()
+    set(CLANG_TIDY_COMMAND "clang-tidy")
+endif()
+
+add_custom_target(
+    clang-tidy
+    COMMAND find \( -path "./src/*" -a -not -path "./src/python/*" -a \( -name "*.cpp"   -not -name "*.test.cpp"\)  \)    -not -name "CircularFifo.hpp" -not -name "ProducerConsumerQueue.hpp" -not -name "VariableSizeClusterFinder.hpp" | xargs -I {} -n 1 -P 10 bash -c "${CLANG_TIDY_COMMAND} --config-file=.clang-tidy -p build {}"
+    WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
+    COMMENT "linting with clang-tidy"
+    VERBATIM
+)
+
+if(AARE_MASTER_PROJECT)
+    set(CMAKE_INSTALL_DIR "share/cmake/${PROJECT_NAME}")
+    set(PROJECT_LIBRARIES slsSupportShared slsDetectorShared slsReceiverShared)
+    include(cmake/package_config.cmake)
+endif()

README.md

@@ -1,2 +1,20 @@
 # aare
 Data analysis library for PSI hybrid detectors
+
+
+## Status
+
+- [ ] Build with CMake on RH8
+- [ ] conda package
+
+
+## Project structure 
+
+include/aare - public headers
+
+
+## Open questions
+
+- How many sub libraries? 
+- Where to place test data? This data is also needed for github actions...
+- What to return to numpy? Our NDArray or a numpy ndarray? Lifetime? 

cmake/FindSphinx.cmake

@@ -0,0 +1,11 @@
+#Look for an executable called sphinx-build
+find_program(SPHINX_EXECUTABLE
+             NAMES sphinx-build sphinx-build-3.6
+             DOC "Path to sphinx-build executable")
+ 
+include(FindPackageHandleStandardArgs)
+ 
+#Handle standard arguments to find_package like REQUIRED and QUIET
+find_package_handle_standard_args(Sphinx
+                                  "Failed to find sphinx-build executable"
+                                  SPHINX_EXECUTABLE)

cmake/helpers.cmake

@@ -0,0 +1,6 @@
+function(default_build_type val)
+if (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
+  message(STATUS "No build type selected, default to Release")
+  set(CMAKE_BUILD_TYPE ${val} CACHE STRING "Build type (default ${val})" FORCE)
+endif()
+endfunction()

cmake/package_config.cmake

@@ -0,0 +1,35 @@
+# This cmake code creates the configuration that is found and used by
+# find_package() of another cmake project
+
+# get lower and upper case project name for the configuration files
+
+# configure and install the configuration files
+include(CMakePackageConfigHelpers)
+
+configure_package_config_file(
+  "${CMAKE_SOURCE_DIR}/cmake/project-config.cmake.in"
+  "${PROJECT_BINARY_DIR}/${PROJECT_NAME_LOWER}-config.cmake"
+  INSTALL_DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/${PROJECT_NAME_LOWER}
+  PATH_VARS CMAKE_INSTALL_DIR)
+
+write_basic_package_version_file(
+  "${PROJECT_BINARY_DIR}/${PROJECT_NAME_LOWER}-config-version.cmake"
+  VERSION ${PROJECT_VERSION}
+  COMPATIBILITY SameMajorVersion
+)
+
+install(FILES
+  "${PROJECT_BINARY_DIR}/${PROJECT_NAME_LOWER}-config.cmake"
+  "${PROJECT_BINARY_DIR}/${PROJECT_NAME_LOWER}-config-version.cmake"
+  COMPONENT devel
+  DESTINATION ${CMAKE_INSTALL_DIR}
+)
+
+
+if (PROJECT_LIBRARIES OR PROJECT_STATIC_LIBRARIES)
+  install(
+    EXPORT "${TARGETS_EXPORT_NAME}"
+    FILE ${PROJECT_NAME_LOWER}-targets.cmake
+    DESTINATION ${CMAKE_INSTALL_DIR}
+    )
+endif ()

cmake/project-config.cmake.in

@@ -0,0 +1,26 @@
+# Config file for @PROJECT_NAME_LOWER@
+#
+# It defines the following variables:
+#
+# @PROJECT_NAME_UPPER@_INCLUDE_DIRS     - include directory
+# @PROJECT_NAME_UPPER@_LIBRARIES        - all dynamic libraries
+# @PROJECT_NAME_UPPER@_STATIC_LIBRARIES - all static libraries
+
+@PACKAGE_INIT@
+
+include(CMakeFindDependencyMacro)
+
+set(SLS_USE_HDF5 "@SLS_USE_HDF5@")
+
+
+find_dependency(Threads)
+
+# Add optional dependencies here
+if (SLS_USE_HDF5)
+    find_dependency(HDF5)
+endif ()
+
+set_and_check(@PROJECT_NAME_UPPER@_CMAKE_INCLUDE_DIRS "@PACKAGE_CMAKE_INSTALL_DIR@")
+
+include("${CMAKE_CURRENT_LIST_DIR}/@TARGETS_EXPORT_NAME@.cmake")
+check_required_components("@PROJECT_NAME@")

conda-recipe/build.sh

@@ -0,0 +1,20 @@
+
+
+mkdir build
+mkdir install
+cd build
+cmake .. \
+      -DCMAKE_PREFIX_PATH=$CONDA_PREFIX \
+      -DCMAKE_INSTALL_PREFIX=install \
+      -DAARE_SYSTEM_LIBRARIES=ON \
+      -DAARE_TESTS=ON \
+      -DAARE_PYTHON_BINDINGS=ON \
+      -DCMAKE_BUILD_TYPE=Release \
+
+     
+NCORES=$(getconf _NPROCESSORS_ONLN)
+echo "Building using: ${NCORES} cores"
+cmake --build . -- -j${NCORES}
+cmake --build . --target install
+
+# CTEST_OUTPUT_ON_FAILURE=1 ctest -j 1

conda-recipe/copy_lib.sh

@@ -0,0 +1,19 @@
+mkdir -p $PREFIX/lib
+mkdir -p $PREFIX/bin
+mkdir -p $PREFIX/include/aare
+
+
+#Shared and static libraries
+cp build/install/lib/* $PREFIX/lib/
+
+#Binaries
+# cp build/install/bin/sls_detector_acquire $PREFIX/bin/.
+# cp build/install/bin/sls_detector_get $PREFIX/bin/.
+# cp build/install/bin/sls_detector_put $PREFIX/bin/.
+# cp build/install/bin/sls_detector_help $PREFIX/bin/.
+# cp build/install/bin/slsReceiver $PREFIX/bin/.
+# cp build/install/bin/slsMultiReceiver $PREFIX/bin/.
+
+
+cp build/install/include/aare/* $PREFIX/include/aare
+cp -rv build/install/share $PREFIX

conda-recipe/meta.yaml

@@ -0,0 +1,104 @@
+
+
+package:
+  name: aare_software
+  version: {{ environ.get('GIT_DESCRIBE_TAG', '') }}
+
+source:
+  - path: ..
+
+build:
+  number: 0
+  binary_relocation: True
+  rpaths:
+    - lib/
+
+requirements:
+  build:
+    - {{ compiler('c') }}
+    - {{compiler('cxx')}}
+    - cmake
+    # - qt 5.*
+    # - xorg-libx11
+    # - xorg-libice
+    # - xorg-libxext
+    # - xorg-libsm
+    # - xorg-libxau
+    # - xorg-libxrender
+    # - xorg-libxfixes
+    # - {{ cdt('mesa-libgl-devel') }}      # [linux]
+    # - {{ cdt('mesa-libegl-devel') }}     # [linux]
+    # - {{ cdt('mesa-dri-drivers') }}      # [linux]
+    # - {{ cdt('libselinux') }}            # [linux]
+    # - {{ cdt('libxdamage') }}            # [linux]
+    # - {{ cdt('libxxf86vm') }}            # [linux]
+    # - expat
+
+  host:
+    # - libstdcxx-ng
+    # - libgcc-ng
+    # - xorg-libx11
+    # - xorg-libice
+    # - xorg-libxext
+    # - xorg-libsm
+    # - xorg-libxau
+    # - xorg-libxrender
+    # - xorg-libxfixes
+    # - expat
+
+  run:
+    # - libstdcxx-ng
+    # - libgcc-ng 
+
+
+outputs:
+  - name: aarelib
+    script: copy_lib.sh
+
+    requirements:
+      build:
+        - {{ compiler('c') }}
+        - {{compiler('cxx')}}
+        - catch2
+        - zstd
+
+        # - libstdcxx-ng
+        # - libgcc-ng
+
+
+      run: 
+        # - libstdcxx-ng
+        # - libgcc-ng
+
+  # - name: aare
+
+  #   script: build_pylib.sh
+
+  #   requirements:
+  #     build:
+  #       - python
+  #       - {{ compiler('c') }}
+  #       - {{compiler('cxx')}}
+  #       - {{ pin_subpackage('slsdetlib', exact=True) }}
+  #       - setuptools
+  #       - pybind11=2.11
+
+  #     host:
+  #       - python
+  #       - {{ pin_subpackage('slsdetlib', exact=True) }}
+  #       - pybind11=2.11
+
+
+  #     run:
+  #       - libstdcxx-ng
+  #       - libgcc-ng
+  #       - python
+  #       - numpy
+  #       - {{ pin_subpackage('slsdetlib', exact=True) }}
+
+
+  #   test:
+  #     imports:
+  #       - slsdet
+
+  

docs/CMakeLists.txt

@@ -0,0 +1,63 @@
+find_package(Doxygen REQUIRED)
+find_package(Sphinx REQUIRED)
+
+#Doxygen
+set(DOXYGEN_IN ${CMAKE_CURRENT_SOURCE_DIR}/Doxyfile.in)
+set(DOXYGEN_OUT ${CMAKE_CURRENT_BINARY_DIR}/Doxyfile)
+configure_file(${DOXYGEN_IN} ${DOXYGEN_OUT} @ONLY)
+
+#Sphinx
+set(SPHINX_SOURCE ${CMAKE_CURRENT_SOURCE_DIR}/src)
+set(SPHINX_BUILD ${CMAKE_CURRENT_BINARY_DIR})
+
+set(SPHINX_SOURCE_FILES
+  src/index.rst
+  src/NDArray.rst
+  src/NDView.rst
+  src/File.rst
+  src/Frame.rst
+  src/Dtype.rst
+  src/ClusterFinder.rst
+  src/Pedestal.rst
+  src/VarClusterFinder.rst
+  src/pyFile.rst
+)
+
+foreach(filename ${SPHINX_SOURCE_FILES})
+configure_file(${CMAKE_CURRENT_SOURCE_DIR}/${filename} 
+              "${SPHINX_BUILD}/${filename}")
+endforeach(filename ${SPHINX_SOURCE_FILES})
+
+configure_file(
+    "${CMAKE_CURRENT_SOURCE_DIR}/conf.py.in"
+    "${SPHINX_BUILD}/conf.py"
+    @ONLY
+)
+
+
+configure_file(
+    "${CMAKE_CURRENT_SOURCE_DIR}/static/extra.css"
+    "${SPHINX_BUILD}/static/css/extra.css"
+    @ONLY
+)
+
+add_custom_target(
+    docs
+    COMMAND ${DOXYGEN_EXECUTABLE} ${DOXYGEN_OUT}
+    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"
+)
+
+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/conf.py.in

@@ -0,0 +1,66 @@
+# Configuration file for the Sphinx documentation builder.
+#
+# This file only contains a selection of the most common options. For a full
+# list see the documentation:
+# http://www.sphinx-doc.org/en/master/config
+
+# -- Path setup --------------------------------------------------------------
+
+# If extensions (or modules to document with autodoc) are in another directory,
+# add these directories to sys.path here. If the directory is relative to the
+# documentation root, use os.path.abspath to make it absolute, like shown here.
+#
+import os
+import sys
+# sys.path.insert(0, os.path.abspath('.'))
+sys.path.insert(0, os.path.abspath('../bin/'))
+#sys.path.insert(0, '/home/l_frojdh/sls/build/bin')
+#sys.path.insert(0, @CMAKE_CURRENT_BINARY_DIR@)
+print(sys.path)
+
+# -- Project information -----------------------------------------------------
+
+project = 'aare'
+copyright = '2024, CPS Detector Group'
+author = 'CPS Detector Group'
+version = '@PROJECT_VERSION@'
+
+# -- General configuration ---------------------------------------------------
+
+# Add any Sphinx extension module names here, as strings. They can be
+# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
+# ones.
+extensions = ['breathe',
+              'sphinx_rtd_theme',
+              'sphinx.ext.autodoc',
+              'sphinx.ext.napoleon',
+             ]
+
+breathe_default_project = "aare"
+napoleon_use_ivar = True
+
+
+# Add any paths that contain templates here, relative to this directory.
+templates_path = ['_templates']
+
+# List of patterns, relative to source directory, that match files and
+# directories to ignore when looking for source files.
+# This pattern also affects html_static_path and html_extra_path.
+exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store']
+
+
+# -- Options for HTML output -------------------------------------------------
+
+# The theme to use for HTML and HTML Help pages.  See the documentation for
+# a list of builtin themes.
+#
+html_theme = "furo"
+
+# Add any paths that contain custom static files (such as style sheets) here,
+# relative to this directory. They are copied after the builtin static files,
+# so a file named "default.css" will overwrite the builtin "default.css".
+html_static_path = ['static']
+
+
+def setup(app):
+    app.add_css_file('css/extra.css')  # may also be an URL

docs/src/ClusterFinder.rst

@@ -0,0 +1,7 @@
+ClusterFinder
+=============
+
+
+.. doxygenclass:: aare::ClusterFinder
+   :members:
+   :undoc-members:

docs/src/Dtype.rst

@@ -0,0 +1,7 @@
+Dtype
+=============
+
+
+.. doxygenclass:: aare::Dtype
+   :members:
+   :undoc-members:

docs/src/File.rst

@@ -0,0 +1,7 @@
+File
+=============
+
+
+.. doxygenclass:: aare::File
+   :members:
+   :undoc-members:

docs/src/Frame.rst

@@ -0,0 +1,7 @@
+Frame
+=============
+
+
+.. doxygenclass:: aare::Frame
+   :members:
+   :undoc-members:

docs/src/NDArray.rst

@@ -0,0 +1,7 @@
+NDArray
+=============
+
+
+.. doxygenclass:: aare::NDArray
+   :members:
+   :undoc-members:

docs/src/NDView.rst

@@ -0,0 +1,7 @@
+NDView
+=============
+
+
+.. doxygenclass:: aare::NDView
+   :members:
+   :undoc-members:

docs/src/Pedestal.rst

@@ -0,0 +1,7 @@
+Pedestal
+=============
+
+
+.. doxygenclass:: aare::Pedestal
+   :members:
+   :undoc-members:

docs/src/VarClusterFinder.rst

@@ -0,0 +1,7 @@
+VarClusterFinder
+====================
+
+
+.. doxygenclass:: aare::VarClusterFinder
+   :members:
+   :undoc-members:

docs/src/index.rst

@@ -0,0 +1,25 @@
+AARE
+==============================================
+
+.. note :: 
+
+    Hello
+
+.. toctree::
+    :caption: C++ API
+    :maxdepth: 1
+    
+    NDArray
+    NDView
+    Frame
+    File
+    Dtype
+    ClusterFinder
+    Pedestal
+    VarClusterFinder
+
+.. toctree::
+    :caption: Python API
+    :maxdepth: 1
+    
+    pyFile

docs/src/pyFile.rst

@@ -0,0 +1,11 @@
+
+File
+========
+
+.. py:currentmodule:: aare
+
+.. autoclass:: File
+    :members:
+    :undoc-members:
+    :show-inheritance:
+    :inherited-members:

docs/static/extra.css

@@ -0,0 +1,4 @@
+/* override table no-wrap */
+.wy-table-responsive table td, .wy-table-responsive table th {
+    white-space: normal;
+}

include/aare/ClusterFileV2.hpp

@@ -0,0 +1,148 @@
+#pragma once
+#include "aare/core/defs.hpp"
+#include <filesystem>
+#include <string>
+#include <fmt/format.h>
+
+namespace aare {
+struct ClusterHeader {
+    int32_t frame_number;
+    int32_t n_clusters;
+    std::string to_string() const {
+        return "frame_number: " + std::to_string(frame_number) + ", n_clusters: " + std::to_string(n_clusters);
+    }
+};
+
+struct ClusterV2_ {
+    int16_t x;
+    int16_t y;
+    std::array<int32_t, 9> data;
+    std::string to_string(bool detailed = false) const {
+        if (detailed) {
+            std::string data_str = "[";
+            for (auto &d : data) {
+                data_str += std::to_string(d) + ", ";
+            }
+            data_str += "]";
+            return "x: " + std::to_string(x) + ", y: " + std::to_string(y) + ", data: " + data_str;
+        }
+        return "x: " + std::to_string(x) + ", y: " + std::to_string(y);
+    }
+};
+
+struct ClusterV2 {
+    ClusterV2_ cluster;
+    int32_t frame_number;
+    std::string to_string() const {
+        return "frame_number: " + std::to_string(frame_number) + ", " + cluster.to_string();
+    }
+};
+
+/**
+ * @brief
+ * important not: fp always points to the clusters header and does not point to individual clusters
+ *
+ */
+class ClusterFileV2 {
+    std::filesystem::path m_fpath; 
+    std::string m_mode;
+    FILE *fp{nullptr};
+
+    void check_open(){
+        if (!fp)
+            throw std::runtime_error(fmt::format("File: {} not open", m_fpath.string()));
+    }
+
+  public:
+    ClusterFileV2(std::filesystem::path const &fpath, std::string const &mode): m_fpath(fpath), m_mode(mode) {
+        if (m_mode != "r" && m_mode != "w")
+            throw std::invalid_argument("mode must be 'r' or 'w'");
+        if (m_mode == "r" && !std::filesystem::exists(m_fpath))
+            throw std::invalid_argument("File does not exist");
+        if (mode == "r") {
+            fp = fopen(fpath.string().c_str(), "rb");
+        } else if (mode == "w") {
+            if (std::filesystem::exists(fpath)) {
+                fp = fopen(fpath.string().c_str(), "r+b");
+            } else {
+                fp = fopen(fpath.string().c_str(), "wb");
+            }
+        }
+        if (fp == nullptr) {
+            throw std::runtime_error("Failed to open file");
+        }
+    }
+    ~ClusterFileV2() { close(); }
+    std::vector<ClusterV2> read() {
+        check_open();
+
+        ClusterHeader header;
+        fread(&header, sizeof(ClusterHeader), 1, fp); 
+        std::vector<ClusterV2_> clusters_(header.n_clusters);
+        fread(clusters_.data(), sizeof(ClusterV2_), header.n_clusters, fp);
+        std::vector<ClusterV2> clusters;
+        for (auto &c : clusters_) {
+            ClusterV2 cluster;
+            cluster.cluster = std::move(c);
+            cluster.frame_number = header.frame_number;
+            clusters.push_back(cluster);
+        }
+
+        return clusters;
+    }
+    std::vector<std::vector<ClusterV2>> read(int n_frames) {
+        std::vector<std::vector<ClusterV2>> clusters;
+        for (int i = 0; i < n_frames; i++) {
+            clusters.push_back(read());
+        }
+        return clusters;
+    }
+
+    size_t write(std::vector<ClusterV2> const &clusters) {
+        check_open();
+        if (m_mode != "w")
+            throw std::runtime_error("File not opened in write mode");
+        if (clusters.empty())
+            return 0;
+
+        ClusterHeader header;
+        header.frame_number = clusters[0].frame_number;
+        header.n_clusters = clusters.size();
+        fwrite(&header, sizeof(ClusterHeader), 1, fp);
+        for (auto &c : clusters) {
+            fwrite(&c.cluster, sizeof(ClusterV2_), 1, fp);
+        }
+        return clusters.size();
+    }
+
+    size_t write(std::vector<std::vector<ClusterV2>> const &clusters) {
+        check_open();
+        if (m_mode != "w") 
+            throw std::runtime_error("File not opened in write mode");
+
+        size_t n_clusters = 0;
+        for (auto &c : clusters) {
+            n_clusters += write(c);
+        }
+        return n_clusters;
+    }
+
+    int seek_to_begin() { return fseek(fp, 0, SEEK_SET); }
+    int seek_to_end() { return fseek(fp, 0, SEEK_END); }
+
+    int32_t frame_number() {
+        auto pos = ftell(fp);
+        ClusterHeader header;
+        fread(&header, sizeof(ClusterHeader), 1, fp);
+        fseek(fp, pos, SEEK_SET);
+        return header.frame_number;
+    }
+
+    void close() {
+        if (fp) {
+            fclose(fp);
+            fp = nullptr;
+        }
+    }
+};
+} // namespace aare

include/aare/ClusterFinder.hpp

@@ -0,0 +1,216 @@
+#pragma once
+#include "aare/Dtype.hpp"
+#include "aare/NDArray.hpp"
+#include "aare/NDView.hpp"
+#include "aare/Pedestal.hpp"
+#include <cstddef>
+
+namespace aare {
+
+/** enum to define the event types */
+enum eventType {
+    PEDESTAL,            /** pedestal */
+    NEIGHBOUR,           /** neighbour i.e. below threshold, but in the cluster of a photon */
+    PHOTON,              /** photon i.e. above threshold */
+    PHOTON_MAX,          /** maximum of a cluster satisfying the photon conditions */
+    NEGATIVE_PEDESTAL,   /** negative value, will not be accounted for as pedestal in order to
+                             avoid drift of the pedestal towards negative values */
+    UNDEFINED_EVENT = -1 /** undefined */
+};
+
+class ClusterFinder {
+  public:
+    ClusterFinder(int cluster_sizeX, int cluster_sizeY, double nSigma = 5.0, double threshold = 0.0)
+        : m_cluster_sizeX(cluster_sizeX), m_cluster_sizeY(cluster_sizeY), m_threshold(threshold), m_nSigma(nSigma) {
+
+        c2 = sqrt((cluster_sizeY + 1) / 2 * (cluster_sizeX + 1) / 2);
+        c3 = sqrt(cluster_sizeX * cluster_sizeY);
+    };
+
+    template <typename FRAME_TYPE, typename PEDESTAL_TYPE>
+    std::vector<Cluster> find_clusters_without_threshold(NDView<FRAME_TYPE, 2> frame, Pedestal<PEDESTAL_TYPE> &pedestal,
+                                                         bool late_update = false) {
+        struct pedestal_update {
+            int x;
+            int y;
+            FRAME_TYPE value;
+        };
+        std::vector<pedestal_update> pedestal_updates;
+
+        std::vector<Cluster> clusters;
+        std::vector<std::vector<eventType>> eventMask;
+        for (int i = 0; i < frame.shape(0); i++) {
+            eventMask.push_back(std::vector<eventType>(frame.shape(1)));
+        }
+        long double val;
+        long double max;
+
+        for (int iy = 0; iy < frame.shape(0); iy++) {
+            for (int ix = 0; ix < frame.shape(1); ix++) {
+                // initialize max and total
+                max = std::numeric_limits<FRAME_TYPE>::min();
+                long double total = 0;
+                eventMask[iy][ix] = PEDESTAL;
+
+                for (short ir = -(m_cluster_sizeY / 2); ir < (m_cluster_sizeY / 2) + 1; ir++) {
+                    for (short ic = -(m_cluster_sizeX / 2); ic < (m_cluster_sizeX / 2) + 1; ic++) {
+                        if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) {
+                            val = frame(iy + ir, ix + ic) - pedestal.mean(iy + ir, ix + ic);
+                            total += val;
+                            if (val > max) {
+                                max = val;
+                            }
+                        }
+                    }
+                }
+                auto rms = pedestal.standard_deviation(iy, ix);
+
+                if (frame(iy, ix) - pedestal.mean(iy, ix) < -m_nSigma * rms) {
+                    eventMask[iy][ix] = NEGATIVE_PEDESTAL;
+                    continue;
+                } else if (max > m_nSigma * rms) {
+                    eventMask[iy][ix] = PHOTON;
+
+                } else if (total > c3 * m_nSigma * rms) {
+                    eventMask[iy][ix] = PHOTON;
+                } else{
+                    if (late_update) {
+                        pedestal_updates.push_back({ix, iy, frame(iy, ix)});
+                    } else {
+                        pedestal.push(iy, ix, frame(iy, ix));
+                    }
+                    continue;
+                }
+                if (eventMask[iy][ix] == PHOTON && (frame(iy, ix) - pedestal.mean(iy, ix)) >= max) {
+                    eventMask[iy][ix] = PHOTON_MAX;
+                    Cluster cluster(m_cluster_sizeX, m_cluster_sizeY, Dtype(typeid(PEDESTAL_TYPE)));
+                    cluster.x = ix;
+                    cluster.y = iy;
+                    short i = 0;
+
+                    for (short ir = -(m_cluster_sizeY / 2); ir < (m_cluster_sizeY / 2) + 1; ir++) {
+                        for (short ic = -(m_cluster_sizeX / 2); ic < (m_cluster_sizeX / 2) + 1; ic++) {
+                            if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) {
+                                PEDESTAL_TYPE tmp = static_cast<PEDESTAL_TYPE>(frame(iy + ir, ix + ic)) -
+                                                    pedestal.mean(iy + ir, ix + ic);
+                                cluster.set<PEDESTAL_TYPE>(i, tmp);
+                                i++;
+                            }
+                        }
+                    }
+                    clusters.push_back(cluster);
+                }
+            }
+        }
+        if (late_update) {
+            for (auto &update : pedestal_updates) {
+                pedestal.push(update.y, update.x, update.value);
+            }
+        }
+        return clusters;
+    }
+
+    template <typename FRAME_TYPE, typename PEDESTAL_TYPE>
+    std::vector<Cluster> find_clusters_with_threshold(NDView<FRAME_TYPE, 2> frame, Pedestal<PEDESTAL_TYPE> &pedestal) {
+        assert(m_threshold > 0);
+        std::vector<Cluster> clusters;
+        std::vector<std::vector<eventType>> eventMask;
+        for (int i = 0; i < frame.shape(0); i++) {
+            eventMask.push_back(std::vector<eventType>(frame.shape(1)));
+        }
+        double tthr, tthr1, tthr2;
+
+        NDArray<FRAME_TYPE, 2> rest({frame.shape(0), frame.shape(1)});
+        NDArray<int, 2> nph({frame.shape(0), frame.shape(1)});
+        // convert to n photons
+        // nph = (frame-pedestal.mean()+0.5*m_threshold)/m_threshold; // can be optimized with expression templates?
+        for (int iy = 0; iy < frame.shape(0); iy++) {
+            for (int ix = 0; ix < frame.shape(1); ix++) {
+                auto val = frame(iy, ix) - pedestal.mean(iy, ix);
+                nph(iy, ix) = (val + 0.5 * m_threshold) / m_threshold;
+                nph(iy, ix) = nph(iy, ix) < 0 ? 0 : nph(iy, ix);
+                rest(iy, ix) = val - nph(iy, ix) * m_threshold;
+            }
+        }
+        // iterate over frame pixels
+        for (int iy = 0; iy < frame.shape(0); iy++) {
+            for (int ix = 0; ix < frame.shape(1); ix++) {
+                eventMask[iy][ix] = PEDESTAL;
+                // initialize max and total
+                FRAME_TYPE max = std::numeric_limits<FRAME_TYPE>::min();
+                long double total = 0;
+                if (rest(iy, ix) <= 0.25 * m_threshold) {
+                    pedestal.push(iy, ix, frame(iy, ix));
+                    continue;
+                }
+                eventMask[iy][ix] = NEIGHBOUR;
+                // iterate over cluster pixels around the current pixel (ix,iy)
+                for (short ir = -(m_cluster_sizeY / 2); ir < (m_cluster_sizeY / 2) + 1; ir++) {
+                    for (short ic = -(m_cluster_sizeX / 2); ic < (m_cluster_sizeX / 2) + 1; ic++) {
+                        if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) {
+                            auto val = frame(iy + ir, ix + ic) - pedestal.mean(iy + ir, ix + ic);
+                            total += val;
+                            if (val > max) {
+                                max = val;
+                            }
+                        }
+                    }
+                }
+
+                auto rms = pedestal.standard_deviation(iy, ix);
+                if (m_nSigma == 0) {
+                    tthr = m_threshold;
+                    tthr1 = m_threshold;
+                    tthr2 = m_threshold;
+                } else {
+                    tthr = m_nSigma * rms;
+                    tthr1 = m_nSigma * rms * c3;
+                    tthr2 = m_nSigma * rms * c2;
+
+                    if (m_threshold > 2 * tthr)
+                        tthr = m_threshold - tthr;
+                    if (m_threshold > 2 * tthr1)
+                        tthr1 = tthr - tthr1;
+                    if (m_threshold > 2 * tthr2)
+                        tthr2 = tthr - tthr2;
+                }
+                if (total > tthr1 || max > tthr) {
+                    eventMask[iy][ix] = PHOTON;
+                    nph(iy, ix) += 1;
+                    rest(iy, ix) -= m_threshold;
+                } else {
+                    pedestal.push(iy, ix, frame(iy, ix));
+                    continue;
+                }
+                if (eventMask[iy][ix] == PHOTON && frame(iy, ix) - pedestal.mean(iy, ix) >= max) {
+                    eventMask[iy][ix] = PHOTON_MAX;
+                    Cluster cluster(m_cluster_sizeX, m_cluster_sizeY, Dtype(typeid(FRAME_TYPE)));
+                    cluster.x = ix;
+                    cluster.y = iy;
+                    short i = 0;
+                    for (short ir = -(m_cluster_sizeY / 2); ir < (m_cluster_sizeY / 2) + 1; ir++) {
+                        for (short ic = -(m_cluster_sizeX / 2); ic < (m_cluster_sizeX / 2) + 1; ic++) {
+                            if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) {
+                                auto tmp = frame(iy + ir, ix + ic) - pedestal.mean(iy + ir, ix + ic);
+                                cluster.set<FRAME_TYPE>(i, tmp);
+                                i++;
+                            }
+                        }
+                    }
+                    clusters.push_back(cluster);
+                }
+            }
+        }
+        return clusters;
+    }
+
+  protected:
+    int m_cluster_sizeX;
+    int m_cluster_sizeY;
+    double m_threshold;
+    double m_nSigma;
+    double c2;
+    double c3;
+};
+
+} // namespace aare

include/aare/Dtype.hpp

@@ -0,0 +1,83 @@
+#pragma once
+#include <cstdint>
+#include <map>
+#include <string>
+#include <typeinfo>
+
+namespace aare {
+
+// The format descriptor is a single character that specifies the type of the data
+// - python documentation: https://docs.python.org/3/c-api/arg.html#numbers
+// - py::format_descriptor<T>::format() (in pybind11) does not return the same format as
+//   written in python.org documentation.
+// - numpy also doesn't use the same format. and also numpy associates the format
+//   with variable bitdepth types. (e.g. long is int64 on linux64 and int32 on win64)
+//   https://numpy.org/doc/stable/reference/arrays.scalars.html
+//
+// github issue discussing this:
+// https://github.com/pybind/pybind11/issues/1908#issuecomment-658358767
+//
+// [IN LINUX] the difference is for int64 (long) and uint64 (unsigned long). The format
+// descriptor is 'q' and 'Q' respectively and in the documentation it is 'l' and 'k'.
+
+// in practice numpy doesn't seem to care when reading buffer info: the library
+// interprets 'q' or 'l' as int64 and 'Q' or 'L' as uint64.
+// for this reason we decided to use the same format descriptor as pybind to avoid
+// any further discrepancies.
+
+// in the following order:
+// int8, uint8, int16, uint16, int32, uint32, int64, uint64, float, double
+const char DTYPE_FORMAT_DSC[] = {'b', 'B', 'h', 'H', 'i', 'I', 'q', 'Q', 'f', 'd'};
+
+// on linux64 & apple
+const char NUMPY_FORMAT_DSC[] = {'b', 'B', 'h', 'H', 'i', 'I', 'l', 'L', 'f', 'd'};
+/**
+ * @brief enum class to define the endianess of the system
+ */
+enum class endian {
+#ifdef _WIN32
+    little = 0,
+    big = 1,
+    native = little
+#else
+    little = __ORDER_LITTLE_ENDIAN__,
+    big = __ORDER_BIG_ENDIAN__,
+    native = __BYTE_ORDER__
+#endif
+};
+
+/**
+ * @brief class to define the data type of the pixels
+ * @note only native endianess is supported
+ */
+class Dtype {
+  public:
+    enum TypeIndex { INT8, UINT8, INT16, UINT16, INT32, UINT32, INT64, UINT64, FLOAT, DOUBLE, ERROR, NONE };
+
+    uint8_t bitdepth() const;
+    size_t bytes() const;
+    std::string format_descr() const { return std::string(1, DTYPE_FORMAT_DSC[static_cast<int>(m_type)]); }
+    std::string numpy_descr() const { return std::string(1, NUMPY_FORMAT_DSC[static_cast<int>(m_type)]); }
+
+    explicit Dtype(const std::type_info &t);
+    explicit Dtype(std::string_view sv);
+    static Dtype from_bitdepth(uint8_t bitdepth);
+
+    // not explicit to allow conversions form enum to DType
+    Dtype(Dtype::TypeIndex ti); // NOLINT
+
+    bool operator==(const Dtype &other) const noexcept;
+    bool operator!=(const Dtype &other) const noexcept;
+    bool operator==(const std::type_info &t) const;
+    bool operator!=(const std::type_info &t) const;
+
+    // bool operator==(DType::TypeIndex ti) const;
+    // bool operator!=(DType::TypeIndex ti) const;
+    std::string to_string() const;
+    void set_type(Dtype::TypeIndex ti) { m_type = ti; }
+
+  private:
+    TypeIndex m_type{TypeIndex::ERROR};
+};
+
+} // namespace aare

include/aare/File.hpp

@@ -0,0 +1,60 @@
+#pragma once
+#include "aare/FileInterface.hpp"
+
+
+namespace aare {
+
+/**
+ * @brief RAII File class for reading and writing image files in various formats
+ * wrapper on a FileInterface to abstract the underlying file format
+ * @note documentation for each function is in the FileInterface class
+ */
+class File {
+  private:
+    FileInterface *file_impl;
+    bool is_npy;
+
+  public:
+    /**
+     * @brief Construct a new File object
+     * @param fname path to the file
+     * @param mode file mode (r, w, a)
+     * @param cfg file configuration
+     * @throws std::runtime_error if the file cannot be opened
+     * @throws std::invalid_argument if the file mode is not supported
+     *
+     */
+    File(const std::filesystem::path &fname, const std::string &mode, const FileConfig &cfg = {});
+    void write(Frame &frame, sls_detector_header header = {});
+    Frame read();
+    Frame iread(size_t frame_number);
+    std::vector<Frame> read(size_t n_frames);
+    void read_into(std::byte *image_buf);
+    void read_into(std::byte *image_buf, size_t n_frames);
+    size_t frame_number(size_t frame_index);
+    size_t bytes_per_frame();
+    size_t pixels_per_frame();
+    void seek(size_t frame_number);
+    size_t tell() const;
+    size_t total_frames() const;
+    size_t rows() const;
+    size_t cols() const;
+    size_t bitdepth() const;
+    size_t bytes_per_pixel() const;
+    void set_total_frames(size_t total_frames);
+    DetectorType detector_type() const;
+    xy geometry() const;
+
+    /**
+     * @brief Move constructor
+     * @param other File object to move from
+     */
+    File(File &&other) noexcept;
+
+    /**
+     * @brief destructor: will only delete the FileInterface object
+     */
+    ~File();
+};
+
+} // namespace aare

include/aare/FileInterface.hpp

@@ -0,0 +1,196 @@
+#pragma once
+#include "aare/Dtype.hpp"
+#include "aare/Frame.hpp"
+#include "aare/defs.hpp"
+
+#include <filesystem>
+#include <vector>
+
+namespace aare {
+
+/**
+ * @brief FileConfig structure to store the configuration of a file
+ * dtype: data type of the file
+ * rows: number of rows in the file
+ * cols: number of columns in the file
+ * geometry: geometry of the file
+ */
+struct FileConfig {
+    aare::Dtype dtype{typeid(uint16_t)};
+    uint64_t rows{};
+    uint64_t cols{};
+    bool operator==(const FileConfig &other) const {
+        return dtype == other.dtype && rows == other.rows && cols == other.cols && geometry == other.geometry &&
+               detector_type == other.detector_type && max_frames_per_file == other.max_frames_per_file;
+    }
+    bool operator!=(const FileConfig &other) const { return !(*this == other); }
+
+    // rawfile specific
+    std::string version{};
+    xy geometry{1, 1};
+    DetectorType detector_type{DetectorType::Unknown};
+    int max_frames_per_file{};
+    size_t total_frames{};
+    std::string to_string() const {
+        return "{ dtype: " + dtype.to_string() + ", rows: " + std::to_string(rows) + ", cols: " + std::to_string(cols) +
+               ", geometry: " + geometry.to_string() + ", detector_type: " + toString(detector_type) +
+               ", max_frames_per_file: " + std::to_string(max_frames_per_file) +
+               ", total_frames: " + std::to_string(total_frames) + " }";
+    }
+};
+
+/**
+ * @brief FileInterface class to define the interface for file operations
+ * @note parent class for NumpyFile and RawFile
+ * @note all functions are pure virtual and must be implemented by the derived classes
+ */
+class FileInterface {
+  public:
+    /**
+     * @brief write a frame to the file
+     * @param frame frame to write
+     * @return void
+     * @throws std::runtime_error if the function is not implemented
+     */
+    // virtual void write(Frame &frame) = 0;
+
+    /**
+     * @brief write a vector of frames to the file
+     * @param frames vector of frames to write
+     * @return void
+     */
+    // virtual void write(std::vector<Frame> &frames) = 0;
+
+    /**
+     * @brief read one frame from the file at the current position
+     * @return Frame
+     */
+    virtual Frame read() = 0;
+
+    /**
+     * @brief read n_frames from the file at the current position
+     * @param n_frames number of frames to read
+     * @return vector of frames
+     */
+    virtual std::vector<Frame> read(size_t n_frames) = 0; // Is this the right interface?
+
+    /**
+     * @brief read one frame from the file at the current position and store it in the provided buffer
+     * @param image_buf buffer to store the frame
+     * @return void
+     */
+    virtual void read_into(std::byte *image_buf) = 0;
+
+    /**
+     * @brief read n_frames from the file at the current position and store them in the provided buffer
+     * @param image_buf buffer to store the frames
+     * @param n_frames number of frames to read
+     * @return void
+     */
+    virtual void read_into(std::byte *image_buf, size_t n_frames) = 0;
+
+    /**
+     * @brief get the frame number at the given frame index
+     * @param frame_index index of the frame
+     * @return frame number
+     */
+    virtual size_t frame_number(size_t frame_index) = 0;
+
+    /**
+     * @brief get the size of one frame in bytes
+     * @return size of one frame
+     */
+    virtual size_t bytes_per_frame() = 0;
+
+    /**
+     * @brief get the number of pixels in one frame
+     * @return number of pixels in one frame
+     */
+    virtual size_t pixels_per_frame() = 0;
+
+    /**
+     * @brief seek to the given frame number
+     * @param frame_number frame number to seek to
+     * @return void
+     */
+    virtual void seek(size_t frame_number) = 0;
+
+    /**
+     * @brief get the current position of the file pointer
+     * @return current position of the file pointer
+     */
+    virtual size_t tell() = 0;
+
+    /**
+     * @brief get the total number of frames in the file
+     * @return total number of frames in the file
+     */
+    virtual size_t total_frames() const = 0;
+    /**
+     * @brief get the number of rows in the file
+     * @return number of rows in the file
+     */
+    virtual size_t rows() const = 0;
+    /**
+     * @brief get the number of columns in the file
+     * @return number of columns in the file
+     */
+    virtual size_t cols() const = 0;
+    /**
+     * @brief get the bitdepth of the file
+     * @return bitdepth of the file
+     */
+    virtual size_t bitdepth() const = 0;
+
+    /**
+     * @brief read one frame from the file at the given frame number
+     * @param frame_number frame number to read
+     * @return frame
+     */
+    Frame iread(size_t frame_number) {
+        auto old_pos = tell();
+        seek(frame_number);
+        Frame tmp = read();
+        seek(old_pos);
+        return tmp;
+    };
+
+    /**
+     * @brief read n_frames from the file starting at the given frame number
+     * @param frame_number frame number to start reading from
+     * @param n_frames number of frames to read
+     * @return vector of frames
+     */
+    std::vector<Frame> iread(size_t frame_number, size_t n_frames) {
+        auto old_pos = tell();
+        seek(frame_number);
+        std::vector<Frame> tmp = read(n_frames);
+        seek(old_pos);
+        return tmp;
+    }
+    DetectorType detector_type() const { return m_type; }
+
+    // function to query the data type of the file
+    /*virtual DataType dtype = 0; */
+
+    virtual ~FileInterface() = default;
+
+    void set_total_frames(size_t total_frames) { m_total_frames = total_frames; }
+
+  protected:
+    std::string m_mode{};
+    std::filesystem::path m_fname{};
+    std::filesystem::path m_base_path{};
+    std::string m_base_name{}, m_ext{};
+    int m_findex{};
+    size_t m_total_frames{};
+    size_t max_frames_per_file{};
+    std::string version{};
+    DetectorType m_type{DetectorType::Unknown};
+    size_t m_rows{};
+    size_t m_cols{};
+    size_t m_bitdepth{};
+    size_t current_frame{};
+};
+
+} // namespace aare

include/aare/Frame.hpp

@@ -0,0 +1,76 @@
+#pragma once
+#include "aare/Dtype.hpp"
+#include "aare/NDArray.hpp"
+#include "aare/defs.hpp"
+
+#include <cstddef>
+#include <cstdint>
+#include <memory>
+#include <vector>
+
+namespace aare {
+
+/**
+ * @brief Frame class to represent a single frame of data
+ * model class
+ * should be able to work with streams coming from files or network
+ */
+class Frame {
+    uint32_t m_rows;
+    uint32_t m_cols;
+    Dtype m_dtype;
+    std::byte *m_data;
+
+  public:
+    Frame(uint32_t rows, uint32_t cols, Dtype dtype);
+    Frame(const std::byte *bytes, uint32_t rows, uint32_t cols, Dtype dtype);
+    ~Frame() noexcept;
+
+    // disable copy and assignment
+    Frame &operator=(const Frame &other)=delete;
+    Frame(const Frame &other)=delete;
+
+    // enable move
+    Frame &operator=(Frame &&other) noexcept;
+    Frame(Frame &&other) noexcept;
+
+    // explicit copy
+    Frame copy() const;
+
+    uint32_t rows() const;
+    uint32_t cols() const;
+    size_t bitdepth() const;
+    Dtype dtype() const;
+    uint64_t size() const;
+    size_t bytes() const;
+    std::byte *data() const;
+
+    std::byte *get(uint32_t row, uint32_t col);
+
+    // TODO! can we, or even want to remove the template?
+    template <typename T> void set(uint32_t row, uint32_t col, T data) {
+        assert(sizeof(T) == m_dtype.bytes());
+        if (row >= m_rows || col >= m_cols) {
+            throw std::out_of_range("Invalid row or column index");
+        }
+        std::memcpy(m_data + (row * m_cols + col) * m_dtype.bytes(), &data, m_dtype.bytes());
+    }
+    template <typename T> T get_t(uint32_t row, uint32_t col) {
+        assert(sizeof(T) == m_dtype.bytes());
+        if (row >= m_rows || col >= m_cols) {
+            throw std::out_of_range("Invalid row or column index");
+        }
+        T data;
+        std::memcpy(&data, m_data + (row * m_cols + col) * m_dtype.bytes(), m_dtype.bytes());
+        return data;
+    }
+    template <typename T> NDView<T, 2> view() {
+        std::array<int64_t, 2> shape = {static_cast<int64_t>(m_rows), static_cast<int64_t>(m_cols)};
+        T *data = reinterpret_cast<T *>(m_data);
+        return NDView<T, 2>(data, shape);
+    }
+
+    template <typename T> NDArray<T> image() { return NDArray<T>(this->view<T>()); }
+};
+
+} // namespace aare

include/aare/NDArray.hpp

@@ -0,0 +1,380 @@
+#pragma once
+/*
+Container holding image data, or a time series of image data in contigious
+memory.
+
+
+TODO! Add expression templates for operators
+
+*/
+#include "aare/NDView.hpp"
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <fmt/format.h>
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <numeric>
+
+namespace aare {
+
+template <typename T, int64_t Ndim = 2> class NDArray {
+  public:
+    NDArray() : shape_(), strides_(c_strides<Ndim>(shape_)), data_(nullptr){};
+
+    explicit NDArray(std::array<int64_t, Ndim> shape)
+        : shape_(shape), strides_(c_strides<Ndim>(shape_)),
+          size_(std::accumulate(shape_.begin(), shape_.end(), 1, std::multiplies<>())), data_(new T[size_]){};
+
+    NDArray(std::array<int64_t, Ndim> shape, T value) : NDArray(shape) { this->operator=(value); }
+
+    /* When constructing from a NDView we need to copy the data since
+    NDArray expect to own its data, and span is just a view*/
+    explicit NDArray(NDView<T, Ndim> span) : NDArray(span.shape()) {
+        std::copy(span.begin(), span.end(), begin());
+        // fmt::print("NDArray(NDView<T, Ndim> span)\n");
+    }
+
+    // Move constructor
+    NDArray(NDArray &&other) noexcept
+        : shape_(other.shape_), strides_(c_strides<Ndim>(shape_)), size_(other.size_), data_(other.data_) {
+
+        other.reset();
+        // fmt::print("NDArray(NDArray &&other)\n");
+    }
+
+    // Copy constructor
+    NDArray(const NDArray &other)
+        : shape_(other.shape_), strides_(c_strides<Ndim>(shape_)), size_(other.size_), data_(new T[size_]) {
+        std::copy(other.data_, other.data_ + size_, data_);
+        // fmt::print("NDArray(const NDArray &other)\n");
+    }
+
+    ~NDArray() { delete[] data_; }
+
+    auto begin() { return data_; }
+    auto end() { return data_ + size_; }
+
+    using value_type = T;
+
+    NDArray &operator=(NDArray &&other) noexcept; // Move assign
+    NDArray &operator=(const NDArray &other);     // Copy assign
+
+    NDArray operator+(const NDArray &other);
+    NDArray &operator+=(const NDArray &other);
+    NDArray operator-(const NDArray &other);
+    NDArray &operator-=(const NDArray &other);
+    NDArray operator*(const NDArray &other);
+    NDArray &operator*=(const NDArray &other);
+    NDArray operator/(const NDArray &other);
+    // NDArray& operator/=(const NDArray& other);
+    template <typename V> NDArray &operator/=(const NDArray<V, Ndim> &other) {
+        // check shape
+        if (shape_ == other.shape()) {
+            for (uint32_t i = 0; i < size_; ++i) {
+                data_[i] /= other(i);
+            }
+            return *this;
+        }
+        throw(std::runtime_error("Shape of NDArray must match"));
+    }
+
+    NDArray<bool, Ndim> operator>(const NDArray &other);
+
+    bool operator==(const NDArray &other) const;
+    bool operator!=(const NDArray &other) const;
+
+    NDArray &operator=(const T & /*value*/);
+    NDArray &operator+=(const T & /*value*/);
+    NDArray operator+(const T & /*value*/);
+    NDArray &operator-=(const T & /*value*/);
+    NDArray operator-(const T & /*value*/);
+    NDArray &operator*=(const T & /*value*/);
+    NDArray operator*(const T & /*value*/);
+    NDArray &operator/=(const T & /*value*/);
+    NDArray operator/(const T & /*value*/);
+
+    NDArray &operator&=(const T & /*mask*/);
+
+    void sqrt() {
+        for (int i = 0; i < size_; ++i) {
+            data_[i] = std::sqrt(data_[i]);
+        }
+    }
+
+    NDArray &operator++(); // pre inc
+
+    template <typename... Ix> std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) {
+        return data_[element_offset(strides_, index...)];
+    }
+
+    template <typename... Ix> std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) const {
+        return data_[element_offset(strides_, index...)];
+    }
+
+    template <typename... Ix> std::enable_if_t<sizeof...(Ix) == Ndim, T> value(Ix... index) {
+        return data_[element_offset(strides_, index...)];
+    }
+
+    T &operator()(int i) { return data_[i]; }
+    const T &operator()(int i) const { return data_[i]; }
+
+    T *data() { return data_; }
+    std::byte *buffer() { return reinterpret_cast<std::byte *>(data_); }
+    uint64_t size() const { return size_; }
+    size_t total_bytes() const { return size_ * sizeof(T); }
+    std::array<int64_t, Ndim> shape() const noexcept { return shape_; }
+    int64_t shape(int64_t i) const noexcept { return shape_[i]; }
+    std::array<int64_t, Ndim> strides() const noexcept { return strides_; }
+    size_t bitdepth() const noexcept { return sizeof(T) * 8; }
+    std::array<int64_t, Ndim> byte_strides() const noexcept {
+        auto byte_strides = strides_;
+        for (auto &val : byte_strides)
+            val *= sizeof(T);
+        return byte_strides;
+        // return strides_;
+    }
+
+    NDView<T, Ndim> span() const { return NDView<T, Ndim>{data_, shape_}; }
+
+    void Print();
+    void Print_all();
+    void Print_some();
+
+    void reset() {
+        data_ = nullptr;
+        size_ = 0;
+        std::fill(shape_.begin(), shape_.end(), 0);
+        std::fill(strides_.begin(), strides_.end(), 0);
+    }
+
+  private:
+    std::array<int64_t, Ndim> shape_;
+    std::array<int64_t, Ndim> strides_;
+    uint64_t size_{};
+    T *data_;
+};
+
+// Move assign
+template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator=(NDArray<T, Ndim> &&other) noexcept {
+    if (this != &other) {
+        delete[] data_;
+        data_ = other.data_;
+        shape_ = other.shape_;
+        size_ = other.size_;
+        strides_ = other.strides_;
+        other.reset();
+    }
+    return *this;
+}
+
+template <typename T, int64_t Ndim> NDArray<T, Ndim> NDArray<T, Ndim>::operator+(const NDArray &other) {
+    NDArray result(*this);
+    result += other;
+    return result;
+}
+template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const NDArray<T, Ndim> &other) {
+    // check shape
+    if (shape_ == other.shape_) {
+        for (uint32_t i = 0; i < size_; ++i) {
+            data_[i] += other.data_[i];
+        }
+        return *this;
+    }
+    throw(std::runtime_error("Shape of ImageDatas must match"));
+}
+
+template <typename T, int64_t Ndim> NDArray<T, Ndim> NDArray<T, Ndim>::operator-(const NDArray &other) {
+    NDArray result{*this};
+    result -= other;
+    return result;
+}
+
+template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator-=(const NDArray<T, Ndim> &other) {
+    // check shape
+    if (shape_ == other.shape_) {
+        for (uint32_t i = 0; i < size_; ++i) {
+            data_[i] -= other.data_[i];
+        }
+        return *this;
+    }
+    throw(std::runtime_error("Shape of ImageDatas must match"));
+}
+template <typename T, int64_t Ndim> NDArray<T, Ndim> NDArray<T, Ndim>::operator*(const NDArray &other) {
+    NDArray result = *this;
+    result *= other;
+    return result;
+}
+
+template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator*=(const NDArray<T, Ndim> &other) {
+    // check shape
+    if (shape_ == other.shape_) {
+        for (uint32_t i = 0; i < size_; ++i) {
+            data_[i] *= other.data_[i];
+        }
+        return *this;
+    }
+    throw(std::runtime_error("Shape of ImageDatas must match"));
+}
+
+template <typename T, int64_t Ndim> NDArray<T, Ndim> NDArray<T, Ndim>::operator/(const NDArray &other) {
+    NDArray result = *this;
+    result /= other;
+    return result;
+}
+
+template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator&=(const T &mask) {
+    for (auto it = begin(); it != end(); ++it)
+        *it &= mask;
+    return *this;
+}
+
+// template <typename T, int64_t Ndim>
+// NDArray<T, Ndim>& NDArray<T, Ndim>::operator/=(const NDArray<T, Ndim>&
+// other)
+// {
+//     //check shape
+//     if (shape_ == other.shape_) {
+//         for (int i = 0; i < size_; ++i) {
+//             data_[i] /= other.data_[i];
+//         }
+//         return *this;
+//     } else {
+//         throw(std::runtime_error("Shape of ImageDatas must match"));
+//     }
+// }
+
+template <typename T, int64_t Ndim> NDArray<bool, Ndim> NDArray<T, Ndim>::operator>(const NDArray &other) {
+    if (shape_ == other.shape_) {
+        NDArray<bool> result{shape_};
+        for (int i = 0; i < size_; ++i) {
+            result(i) = (data_[i] > other.data_[i]);
+        }
+        return result;
+    }
+    throw(std::runtime_error("Shape of ImageDatas must match"));
+}
+
+template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator=(const NDArray<T, Ndim> &other) {
+    if (this != &other) {
+        delete[] data_;
+        shape_ = other.shape_;
+        strides_ = other.strides_;
+        size_ = other.size_;
+        data_ = new T[size_];
+        std::copy(other.data_, other.data_ + size_, data_);
+    }
+    return *this;
+}
+
+template <typename T, int64_t Ndim> bool NDArray<T, Ndim>::operator==(const NDArray<T, Ndim> &other) const {
+    if (shape_ != other.shape_)
+        return false;
+
+    for (uint32_t i = 0; i != size_; ++i)
+        if (data_[i] != other.data_[i])
+            return false;
+
+    return true;
+}
+
+template <typename T, int64_t Ndim> bool NDArray<T, Ndim>::operator!=(const NDArray<T, Ndim> &other) const {
+    return !((*this) == other);
+}
+template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator++() {
+    for (uint32_t i = 0; i < size_; ++i)
+        data_[i] += 1;
+    return *this;
+}
+template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator=(const T &value) {
+    std::fill_n(data_, size_, value);
+    return *this;
+}
+
+template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const T &value) {
+    for (uint32_t i = 0; i < size_; ++i)
+        data_[i] += value;
+    return *this;
+}
+
+template <typename T, int64_t Ndim> NDArray<T, Ndim> NDArray<T, Ndim>::operator+(const T &value) {
+    NDArray result = *this;
+    result += value;
+    return result;
+}
+template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator-=(const T &value) {
+    for (uint32_t i = 0; i < size_; ++i)
+        data_[i] -= value;
+    return *this;
+}
+template <typename T, int64_t Ndim> NDArray<T, Ndim> NDArray<T, Ndim>::operator-(const T &value) {
+    NDArray result = *this;
+    result -= value;
+    return result;
+}
+
+template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator/=(const T &value) {
+    for (uint32_t i = 0; i < size_; ++i)
+        data_[i] /= value;
+    return *this;
+}
+template <typename T, int64_t Ndim> NDArray<T, Ndim> NDArray<T, Ndim>::operator/(const T &value) {
+    NDArray result = *this;
+    result /= value;
+    return result;
+}
+template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator*=(const T &value) {
+    for (uint32_t i = 0; i < size_; ++i)
+        data_[i] *= value;
+    return *this;
+}
+template <typename T, int64_t Ndim> NDArray<T, Ndim> NDArray<T, Ndim>::operator*(const T &value) {
+    NDArray result = *this;
+    result *= value;
+    return result;
+}
+template <typename T, int64_t Ndim> void NDArray<T, Ndim>::Print() {
+    if (shape_[0] < 20 && shape_[1] < 20)
+        Print_all();
+    else
+        Print_some();
+}
+template <typename T, int64_t Ndim> void NDArray<T, Ndim>::Print_all() {
+    for (auto row = 0; row < shape_[0]; ++row) {
+        for (auto col = 0; col < shape_[1]; ++col) {
+            std::cout << std::setw(3);
+            std::cout << (*this)(row, col) << " ";
+        }
+        std::cout << "\n";
+    }
+}
+template <typename T, int64_t Ndim> void NDArray<T, Ndim>::Print_some() {
+    for (auto row = 0; row < 5; ++row) {
+        for (auto col = 0; col < 5; ++col) {
+            std::cout << std::setw(7);
+            std::cout << (*this)(row, col) << " ";
+        }
+        std::cout << "\n";
+    }
+}
+
+template <typename T, int64_t Ndim> void save(NDArray<T, Ndim> &img, std::string &pathname) {
+    std::ofstream f;
+    f.open(pathname, std::ios::binary);
+    f.write(img.buffer(), img.size() * sizeof(T));
+    f.close();
+}
+
+template <typename T, int64_t Ndim>
+NDArray<T, Ndim> load(const std::string &pathname, std::array<int64_t, Ndim> shape) {
+    NDArray<T, Ndim> img{shape};
+    std::ifstream f;
+    f.open(pathname, std::ios::binary);
+    f.read(img.buffer(), img.size() * sizeof(T));
+    f.close();
+    return img;
+}
+
+} // namespace aare

include/aare/NDView.hpp

@@ -0,0 +1,159 @@
+#pragma once
+#include <algorithm>
+#include <array>
+#include <cassert>
+#include <cstdint>
+#include <functional>
+#include <iomanip>
+#include <iostream>
+#include <numeric>
+#include <stdexcept>
+#include <vector>
+namespace aare {
+
+template <int64_t Ndim> using Shape = std::array<int64_t, Ndim>;
+
+// TODO! fix mismatch between signed and unsigned
+template <int64_t Ndim> Shape<Ndim> make_shape(const std::vector<size_t> &shape) {
+    if (shape.size() != Ndim)
+        throw std::runtime_error("Shape size mismatch");
+    Shape<Ndim> arr;
+    std::copy_n(shape.begin(), Ndim, arr.begin());
+    return arr;
+}
+
+template <int64_t Dim = 0, typename Strides> int64_t element_offset(const Strides & /*unused*/) { return 0; }
+
+template <int64_t Dim = 0, typename Strides, typename... Ix>
+int64_t element_offset(const Strides &strides, int64_t i, Ix... index) {
+    return i * strides[Dim] + element_offset<Dim + 1>(strides, index...);
+}
+
+template <int64_t Ndim> std::array<int64_t, Ndim> c_strides(const std::array<int64_t, Ndim> &shape) {
+    std::array<int64_t, Ndim> strides{};
+    std::fill(strides.begin(), strides.end(), 1);
+    for (int64_t i = Ndim - 1; i > 0; --i) {
+        strides[i - 1] = strides[i] * shape[i];
+    }
+    return strides;
+}
+
+template <int64_t Ndim> std::array<int64_t, Ndim> make_array(const std::vector<int64_t> &vec) {
+    assert(vec.size() == Ndim);
+    std::array<int64_t, Ndim> arr{};
+    std::copy_n(vec.begin(), Ndim, arr.begin());
+    return arr;
+}
+
+template <typename T, int64_t Ndim = 2> class NDView {
+  public:
+    NDView() = default;
+    ~NDView() = default;
+    NDView(const NDView &) = default;
+    NDView(NDView &&) = default;
+
+    NDView(T *buffer, std::array<int64_t, Ndim> shape)
+        : buffer_(buffer), strides_(c_strides<Ndim>(shape)), shape_(shape),
+          size_(std::accumulate(std::begin(shape), std::end(shape), 1, std::multiplies<>())) {}
+
+    // NDView(T *buffer, const std::vector<int64_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 {
+        return buffer_[element_offset(strides_, index...)];
+    }
+
+    uint64_t size() const { return size_; }
+    size_t total_bytes() const { return size_ * sizeof(T); }
+    std::array<int64_t, Ndim> strides() const noexcept { return strides_; }
+
+    T *begin() { return buffer_; }
+    T *end() { return buffer_ + size_; }
+    T &operator()(int64_t i) const { return buffer_[i]; }
+    T &operator[](int64_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 (buffer_[i] != other.buffer_[i])
+                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) { return elemenwise(val, std::multiplies<T>()); }
+    NDView &operator/=(const T val) { return elemenwise(val, std::divides<T>()); }
+
+    NDView &operator/=(const NDView &other) { return elemenwise(other, std::divides<T>()); }
+
+    NDView &operator=(const T val) {
+        for (auto it = begin(); it != end(); ++it)
+            *it = val;
+        return *this;
+    }
+
+    NDView &operator=(const NDView &other) {
+        if (this == &other)
+            return *this;
+        shape_ = other.shape_;
+        strides_ = other.strides_;
+        size_ = other.size_;
+        buffer_ = other.buffer_;
+        return *this;
+    }
+
+    NDView &operator=(NDView &&other) noexcept {
+        if (this == &other)
+            return *this;
+        shape_ = std::move(other.shape_);
+        strides_ = std::move(other.strides_);
+        size_ = other.size_;
+        buffer_ = other.buffer_;
+        other.buffer_ = nullptr;
+        return *this;
+    }
+
+    auto &shape() { return shape_; }
+    auto shape(int64_t i) const { return shape_[i]; }
+
+    T *data() { return buffer_; }
+    void print_all() const;
+
+  private:
+    T *buffer_{nullptr};
+    std::array<int64_t, Ndim> strides_{};
+    std::array<int64_t, Ndim> shape_{};
+    uint64_t size_{};
+
+    template <class BinaryOperation> NDView &elemenwise(T val, BinaryOperation op) {
+        for (uint64_t i = 0; i != size_; ++i) {
+            buffer_[i] = op(buffer_[i], val);
+        }
+        return *this;
+    }
+    template <class BinaryOperation> NDView &elemenwise(const NDView &other, BinaryOperation op) {
+        for (uint64_t i = 0; i != size_; ++i) {
+            buffer_[i] = op(buffer_[i], other.buffer_[i]);
+        }
+        return *this;
+    }
+};
+template <typename T, int64_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) {
+            std::cout << std::setw(3);
+            std::cout << (*this)(row, col) << " ";
+        }
+        std::cout << "\n";
+    }
+}
+
+} // namespace aare

include/aare/NumpyFile.hpp

@@ -0,0 +1,112 @@
+#pragma once
+#include "aare/Dtype.hpp"
+#include "aare/defs.hpp"
+#include "aare/FileInterface.hpp"
+#include "aare/NumpyHelpers.hpp"
+
+
+#include <filesystem>
+#include <iostream>
+#include <numeric>
+
+namespace aare {
+
+
+
+/**
+ * @brief NumpyFile class to read and write numpy files
+ * @note derived from FileInterface
+ * @note implements all the pure virtual functions from FileInterface
+ * @note documentation for the functions can also be found in the FileInterface class
+ */
+class NumpyFile : public FileInterface {
+
+  public:
+    /**
+     * @brief NumpyFile constructor
+     * @param fname path to the numpy file
+     * @param mode file mode (r, w)
+     * @param cfg file configuration
+     */
+    explicit NumpyFile(const std::filesystem::path &fname, const std::string &mode = "r", FileConfig cfg = {});
+
+    void write(Frame &frame);
+    Frame read() override { return get_frame(this->current_frame++); }
+
+    std::vector<Frame> read(size_t n_frames) override;
+    void read_into(std::byte *image_buf) override { return get_frame_into(this->current_frame++, image_buf); }
+    void read_into(std::byte *image_buf, size_t n_frames) override;
+    size_t frame_number(size_t frame_index) override { return frame_index; };
+    size_t bytes_per_frame() override;
+    size_t pixels_per_frame() override;
+    void seek(size_t frame_number) override { this->current_frame = frame_number; }
+    size_t tell() override { return this->current_frame; }
+    size_t total_frames() const override { return m_header.shape[0]; }
+    size_t rows() const override { return m_header.shape[1]; }
+    size_t cols() const override { return m_header.shape[2]; }
+    size_t bitdepth() const override { return m_header.dtype.bitdepth(); }
+
+    /**
+     * @brief get the data type of the numpy file
+     * @return DType
+     */
+    Dtype dtype() const { return m_header.dtype; }
+
+    /**
+     * @brief get the shape of the numpy file
+     * @return vector of type size_t
+     */
+    std::vector<size_t> shape() const { return m_header.shape; }
+
+    /**
+     * @brief load the numpy file into an NDArray
+     * @tparam T data type of the NDArray
+     * @tparam NDim number of dimensions of the NDArray
+     * @return NDArray<T, NDim>
+     */
+    template <typename T, size_t NDim> NDArray<T, NDim> load() {
+        NDArray<T, NDim> arr(make_shape<NDim>(m_header.shape));
+        if (fseek(fp, static_cast<int64_t>(header_size), SEEK_SET)) {
+            throw std::runtime_error(LOCATION + "Error seeking to the start of the data");
+        }
+        size_t rc = fread(arr.data(), sizeof(T), arr.size(), fp);
+        if (rc != static_cast<size_t>(arr.size())) {
+            throw std::runtime_error(LOCATION + "Error reading data from file");
+        }
+        return arr;
+    }
+    template <typename A, typename TYPENAME, A Ndim> void write(NDView<TYPENAME, Ndim> &frame) {
+        write_impl(frame.data(), frame.total_bytes());
+    }
+    template <typename A, typename TYPENAME, A Ndim> void write(NDArray<TYPENAME, Ndim> &frame) {
+        write_impl(frame.data(), frame.total_bytes());
+    }
+    template <typename A, typename TYPENAME, A Ndim> void write(NDView<TYPENAME, Ndim> &&frame) {
+        write_impl(frame.data(), frame.total_bytes());
+    }
+    template <typename A, typename TYPENAME, A Ndim> void write(NDArray<TYPENAME, Ndim> &&frame) {
+        write_impl(frame.data(), frame.total_bytes());
+    }
+
+    ~NumpyFile() noexcept override;
+
+  private:
+    FILE *fp = nullptr;
+    size_t initial_header_len = 0;
+    size_t current_frame{};
+    uint32_t header_len{};
+    uint8_t header_len_size{};
+    size_t header_size{};
+    NumpyHeader m_header;
+    uint8_t major_ver_{};
+    uint8_t minor_ver_{};
+    size_t m_bytes_per_frame{};
+    size_t m_pixels_per_frame{};
+
+    void load_metadata();
+    void get_frame_into(size_t /*frame_number*/, std::byte * /*image_buf*/);
+    Frame get_frame(size_t frame_number);
+    void write_impl(void *data, uint64_t size);
+};
+
+} // namespace aare

include/aare/NumpyHelpers.hpp

@@ -0,0 +1,55 @@
+
+#pragma once
+#include <algorithm>
+#include <array>
+#include <filesystem>
+#include <fstream>
+#include <iostream>
+#include <numeric>
+#include <sstream>
+#include <string>
+#include <unordered_map>
+#include <vector>
+
+#include "aare/Dtype.hpp"
+#include "aare/defs.hpp"
+
+namespace aare {
+
+struct NumpyHeader {
+    Dtype dtype{aare::Dtype::ERROR};
+    bool fortran_order{false};
+    std::vector<size_t> shape{};
+    std::string to_string() const;
+};
+
+namespace NumpyHelpers {
+
+const constexpr std::array<char, 6> magic_str{'\x93', 'N', 'U', 'M', 'P', 'Y'};
+const uint8_t magic_string_length{6};
+
+std::string parse_str(const std::string &in);
+/**
+  Removes leading and trailing whitespaces
+  */
+std::string trim(const std::string &str);
+
+std::vector<std::string> parse_tuple(std::string in);
+
+bool parse_bool(const std::string &in);
+
+std::string get_value_from_map(const std::string &mapstr);
+
+std::unordered_map<std::string, std::string> parse_dict(std::string in, const std::vector<std::string> &keys);
+
+template <typename T, size_t N> bool in_array(T val, const std::array<T, N> &arr) {
+    return std::find(std::begin(arr), std::end(arr), val) != std::end(arr);
+}
+bool is_digits(const std::string &str);
+
+aare::Dtype parse_descr(std::string typestring);
+size_t write_header(const std::filesystem::path &fname, const NumpyHeader &header);
+size_t write_header(std::ostream &out, const NumpyHeader &header);
+
+} // namespace NumpyHelpers
+} // namespace aare

include/aare/Pedestal.hpp

@@ -0,0 +1,121 @@
+#pragma once
+#include "aare/Frame.hpp"
+#include "aare/NDArray.hpp"
+#include "aare/NDView.hpp"
+#include <cstddef>
+
+namespace aare {
+
+template <typename SUM_TYPE = double> class Pedestal {
+  public:
+    Pedestal(uint32_t rows, uint32_t cols, uint32_t n_samples = 1000)
+        : m_rows(rows), m_cols(cols), m_freeze(false), m_samples(n_samples),           m_cur_samples(NDArray<uint32_t, 2>({rows, cols}, 0)),m_sum(NDArray<SUM_TYPE, 2>({rows, cols})),
+ m_sum2(NDArray<SUM_TYPE, 2>({rows, cols})) {
+        assert(rows > 0 && cols > 0 && n_samples > 0);
+        m_sum = 0;
+        m_sum2 = 0;
+    }
+    ~Pedestal() = default;
+
+    NDArray<SUM_TYPE, 2> mean() {
+        NDArray<SUM_TYPE, 2> mean_array({m_rows, m_cols});
+        for (uint32_t i = 0; i < m_rows * m_cols; i++) {
+            mean_array(i / m_cols, i % m_cols) = mean(i / m_cols, i % m_cols);
+        }
+        return mean_array;
+    }
+
+    NDArray<SUM_TYPE, 2> variance() {
+        NDArray<SUM_TYPE, 2> variance_array({m_rows, m_cols});
+        for (uint32_t i = 0; i < m_rows * m_cols; i++) {
+            variance_array(i / m_cols, i % m_cols) = variance(i / m_cols, i % m_cols);
+        }
+        return variance_array;
+    }
+
+    NDArray<SUM_TYPE, 2> standard_deviation() {
+        NDArray<SUM_TYPE, 2> standard_deviation_array({m_rows, m_cols});
+        for (uint32_t i = 0; i < m_rows * m_cols; i++) {
+            standard_deviation_array(i / m_cols, i % m_cols) = standard_deviation(i / m_cols, i % m_cols);
+        }
+
+        return standard_deviation_array;
+    }
+    void clear() {
+        for (uint32_t i = 0; i < m_rows * m_cols; i++) {
+            clear(i / m_cols, i % m_cols);
+        }
+    }
+
+    /*
+     * index level operations
+     */
+    SUM_TYPE mean(const uint32_t row, const uint32_t col) const {
+        if (m_cur_samples(row, col) == 0) {
+            return 0.0;
+        }
+        return m_sum(row, col) / m_cur_samples(row, col);
+    }
+    SUM_TYPE variance(const uint32_t row, const uint32_t col) const {
+        if (m_cur_samples(row, col) == 0) {
+            return 0.0;
+        }
+        return m_sum2(row, col) / m_cur_samples(row, col) - mean(row, col) * mean(row, col);
+    }
+    SUM_TYPE standard_deviation(const uint32_t row, const uint32_t col) const { return std::sqrt(variance(row, col)); }
+
+    void clear(const uint32_t row, const uint32_t col) {
+        m_sum(row, col) = 0;
+        m_sum2(row, col) = 0;
+        m_cur_samples(row, col) = 0;
+    }
+    // frame level operations
+    template <typename T> void push(NDView<T, 2> frame) {
+        assert(frame.size() == m_rows * m_cols);
+        // TODO: test the effect of #pragma omp parallel for
+        for (uint32_t index = 0; index < m_rows * m_cols; index++) {
+            push<T>(index / m_cols, index % m_cols, frame(index));
+        }
+    }
+    template <typename T> void push(Frame &frame) {
+        assert(frame.rows() == static_cast<size_t>(m_rows) && frame.cols() == static_cast<size_t>(m_cols));
+        push<T>(frame.view<T>());
+    }
+
+    // getter functions
+    inline uint32_t rows() const { return m_rows; }
+    inline uint32_t cols() const { return m_cols; }
+    inline uint32_t n_samples() const { return m_samples; }
+    inline NDArray<uint32_t, 2> cur_samples() const { return m_cur_samples; }
+    inline NDArray<SUM_TYPE, 2> get_sum() const { return m_sum; }
+    inline NDArray<SUM_TYPE, 2> get_sum2() const { return m_sum2; }
+
+    // pixel level operations (should be refactored to allow users to implement their own pixel level operations)
+    template <typename T> inline void push(const uint32_t row, const uint32_t col, const T val_) {
+        if (m_freeze) {
+            return;
+        }
+        SUM_TYPE val = static_cast<SUM_TYPE>(val_);
+        const uint32_t idx = index(row, col);
+        if (m_cur_samples(idx) < m_samples) {
+            m_sum(idx) += val;
+            m_sum2(idx) += val * val;
+            m_cur_samples(idx)++;
+        } else {
+            m_sum(idx) += val - m_sum(idx) / m_cur_samples(idx);
+            m_sum2(idx) += val * val - m_sum2(idx) / m_cur_samples(idx);
+        }
+    }
+    inline uint32_t index(const uint32_t row, const uint32_t col) const { return row * m_cols + col; };
+    void set_freeze(bool freeze) { m_freeze = freeze; }
+
+  private:
+    uint32_t m_rows;
+    uint32_t m_cols;
+    bool m_freeze;
+    uint32_t m_samples;
+    NDArray<uint32_t, 2> m_cur_samples;
+    NDArray<SUM_TYPE, 2> m_sum;
+    NDArray<SUM_TYPE, 2> m_sum2;
+};
+} // namespace aare

include/aare/RawFile.hpp

@@ -0,0 +1,186 @@
+#pragma once
+#include "aare/Frame.hpp"
+#include "aare/FileInterface.hpp"
+#include "aare/SubFile.hpp"
+
+namespace aare {
+
+struct ModuleConfig {
+    int module_gap_row{};
+    int module_gap_col{};
+
+    bool operator==(const ModuleConfig &other) const {
+        if (module_gap_col != other.module_gap_col)
+            return false;
+        if (module_gap_row != other.module_gap_row)
+            return false;
+        return true;
+    }
+};
+
+/**
+ * @brief RawFile class to read .raw and .json files
+ * @note derived from FileInterface
+ * @note documentation can also be found in the FileInterface class
+ */
+class RawFile : public FileInterface {
+  public:
+    /**
+     * @brief RawFile constructor
+     * @param fname path to the file
+     * @param mode file mode (r, w)
+     * @param cfg file configuration
+     */
+    explicit RawFile(const std::filesystem::path &fname, const std::string &mode = "r",
+                     const FileConfig &config = FileConfig{});
+
+    /**
+     * @brief write function is not implemented for RawFile
+     * @param frame frame to write
+     */
+    void write(Frame &frame, sls_detector_header header);
+    Frame read() override { return get_frame(this->current_frame++); };
+    std::vector<Frame> read(size_t n_frames) override;
+    void read_into(std::byte *image_buf) override { return get_frame_into(this->current_frame++, image_buf); };
+    void read_into(std::byte *image_buf, size_t n_frames) override;
+    size_t frame_number(size_t frame_index) override;
+
+    /**
+     * @brief get the number of bytess per frame
+     * @return size of one frame in bytes
+     */
+    size_t bytes_per_frame() override { return m_rows * m_cols * m_bitdepth / 8; }
+
+    /**
+     * @brief get the number of pixels in the frame
+     * @return number of pixels
+     */
+    size_t pixels_per_frame() override { return m_rows * m_cols; }
+
+    // goto frame index
+    void seek(size_t frame_index) override {
+        // check if the frame number is greater than the total frames
+        // if frame_number == total_frames, then the next read will throw an error
+        if (frame_index > this->total_frames()) {
+            throw std::runtime_error(
+                fmt::format("frame number {} is greater than total frames {}", frame_index, m_total_frames));
+        }
+        this->current_frame = frame_index;
+    };
+
+    // return the position of the file pointer (in number of frames)
+    size_t tell() override { return this->current_frame; };
+
+    /**
+     * @brief check if the file is a master file
+     * @param fpath path to the file
+     */
+    static bool is_master_file(const std::filesystem::path &fpath);
+
+    /**
+     * @brief set the module gap row and column
+     * @param row gap between rows
+     * @param col gap between columns
+     */
+    inline void set_config(int row, int col) {
+        cfg.module_gap_row = row;
+        cfg.module_gap_col = col;
+    }
+    // TODO! Deal with fast quad and missing files
+
+    /**
+     * @brief get the number of subfiles for the RawFile
+     * @return number of subfiles
+     */
+    void find_number_of_subfiles();
+
+    /**
+     * @brief get the master file name path for the RawFile
+     * @return path to the master file
+     */
+    inline std::filesystem::path master_fname();
+    /**
+     * @brief get the data file name path for the RawFile with the given module id and file id
+     * @param mod_id module id
+     * @param file_id file id
+     * @return path to the data file
+     */
+    inline std::filesystem::path data_fname(size_t mod_id, size_t file_id);
+
+    /**
+     * @brief destructor: will delete the subfiles
+     */
+    ~RawFile() noexcept override;
+
+    size_t total_frames() const override { return m_total_frames; }
+    size_t rows() const override { return m_rows; }
+    size_t cols() const override { return m_cols; }
+    size_t bitdepth() const override { return m_bitdepth; }
+    xy geometry() { return m_geometry; }
+
+  private:
+    void write_master_file();
+    /**
+     * @brief read the frame at the given frame index into the image buffer
+     * @param frame_number frame number to read
+     * @param image_buf buffer to store the frame
+     */
+    void get_frame_into(size_t frame_index, std::byte *frame_buffer);
+
+    /**
+     * @brief get the frame at the given frame index
+     * @param frame_number frame number to read
+     * @return Frame
+     */
+    Frame get_frame(size_t frame_index);
+
+    /**
+     * @brief parse the file name to get the extension, base name and index
+     */
+    void parse_fname();
+
+    /**
+     * @brief parse the metadata from the file
+     */
+    void parse_metadata();
+
+    /**
+     * @brief parse the metadata of a .raw file
+     */
+    void parse_raw_metadata();
+
+    /**
+     * @brief parse the metadata of a .json file
+     */
+    void parse_json_metadata();
+
+    /**
+     * @brief finds the geometry of the file
+     */
+    void find_geometry();
+
+    /**
+     * @brief read the header of the file
+     * @param fname path to the data subfile
+     * @return sls_detector_header
+     */
+    static sls_detector_header read_header(const std::filesystem::path &fname);
+
+    /**
+     * @brief open the subfiles
+     */
+    void open_subfiles();
+    void parse_config(const FileConfig &config);
+
+    size_t n_subfiles{};
+    size_t n_subfile_parts{};
+    std::vector<std::vector<SubFile *>> subfiles;
+    size_t subfile_rows{}, subfile_cols{};
+    xy m_geometry{};
+    std::vector<xy> positions;
+    ModuleConfig cfg{0, 0};
+    TimingMode timing_mode{};
+    bool quad{false};
+};
+
+} // namespace aare

include/aare/SubFile.hpp

@@ -0,0 +1,77 @@
+#pragma once
+#include "aare/Frame.hpp"
+#include "aare/defs.hpp"
+
+#include <cstdint>
+#include <filesystem>
+#include <map>
+#include <variant>
+
+namespace aare {
+
+/**
+ * @brief Class to read a subfile from a RawFile
+ */
+class SubFile {
+  public:
+    size_t write_part(std::byte *buffer, sls_detector_header header, size_t frame_index);
+    /**
+     * @brief SubFile constructor
+     * @param fname path to the subfile
+     * @param detector detector type
+     * @param rows number of rows in the subfile
+     * @param cols number of columns in the subfile
+     * @param bitdepth bitdepth of the subfile
+     * @throws std::invalid_argument if the detector,type pair is not supported
+     */
+    SubFile(const std::filesystem::path &fname, DetectorType detector, size_t rows, size_t cols, size_t bitdepth,
+            const std::string &mode = "r");
+
+    /**
+     * @brief read the subfile into a buffer
+     * @param buffer pointer to the buffer to read the data into
+     * @return number of bytes read
+     */
+    size_t read_impl_normal(std::byte *buffer);
+
+    /**
+     * @brief read the subfile into a buffer with the bytes flipped
+     * @param buffer pointer to the buffer to read the data into
+     * @return number of bytes read
+     */
+    template <typename DataType> size_t read_impl_flip(std::byte *buffer);
+
+    /**
+     * @brief read the subfile into a buffer with the bytes reordered
+     * @param buffer pointer to the buffer to read the data into
+     * @return number of bytes read
+     */
+    template <typename DataType> size_t read_impl_reorder(std::byte *buffer);
+
+    /**
+     * @brief read the subfile into a buffer with the bytes reordered and flipped
+     * @param buffer pointer to the buffer to read the data into
+     * @param frame_number frame number to read
+     * @return number of bytes read
+     */
+    size_t get_part(std::byte *buffer, size_t frame_index);
+    size_t frame_number(size_t frame_index);
+
+    // TODO: define the inlines as variables and assign them in constructor
+    inline size_t bytes_per_part() const { return (m_bitdepth / 8) * m_rows * m_cols; }
+    inline size_t pixels_per_part() const { return m_rows * m_cols; }
+
+    ~SubFile();
+
+  protected:
+    FILE *fp = nullptr;
+    size_t m_bitdepth;
+    std::filesystem::path m_fname;
+    size_t m_rows{};
+    size_t m_cols{};
+    std::string m_mode;
+    size_t n_frames{};
+    int m_sub_file_index_{};
+};
+
+} // namespace aare

include/aare/VarClusterFinder.hpp

@@ -0,0 +1,307 @@
+#pragma once
+
+#include <algorithm>
+#include <map>
+#include <unordered_map>
+#include <vector>
+
+#include "aare/NDArray.hpp"
+
+const int MAX_CLUSTER_SIZE = 200;
+namespace aare {
+
+template <typename T> class VarClusterFinder {
+  public:
+    struct Hit {
+        int16_t size{};
+        int16_t row{};
+        int16_t col{};
+        uint16_t reserved{}; // for alignment
+        T energy{};
+        T max{};
+
+        // std::vector<int16_t> rows{};
+        // std::vector<int16_t> cols{};
+        int16_t rows[MAX_CLUSTER_SIZE] = {0};
+        int16_t cols[MAX_CLUSTER_SIZE] = {0};
+        double enes[MAX_CLUSTER_SIZE] = {0};
+    };
+
+  private:
+    const std::array<int64_t, 2> shape_;
+    NDView<T, 2> original_;
+    NDArray<int, 2> labeled_;
+    NDArray<int, 2> peripheral_labeled_;
+    NDArray<bool, 2> binary_; // over threshold flag
+    T threshold_;
+    NDView<T, 2> noiseMap;
+    bool use_noise_map = false;
+    int peripheralThresholdFactor_ = 5;
+    int current_label;
+    const std::array<int, 4> di{{0, -1, -1, -1}};              // row ### 8-neighbour by scaning from left to right
+    const std::array<int, 4> dj{{-1, -1, 0, 1}};               // col ### 8-neighbour by scaning from top to bottom
+    const std::array<int, 8> di_{{0, 0, -1, 1, -1, 1, -1, 1}}; // row
+    const std::array<int, 8> dj_{{-1, 1, 0, 0, 1, -1, -1, 1}}; // col
+    std::map<int, int> child;                                  // heirachy: key: child; val: parent
+    std::unordered_map<int, Hit> h_size;
+    std::vector<Hit> hits;
+    // std::vector<std::vector<int16_t>> row
+    int check_neighbours(int i, int j);
+
+  public:
+    VarClusterFinder(image_shape shape, T threshold)
+        : shape_(shape), labeled_(shape, 0), peripheral_labeled_(shape, 0), binary_(shape), threshold_(threshold) {
+        hits.reserve(2000);
+    }
+
+    NDArray<int, 2> labeled() { return labeled_; }
+
+    void set_noiseMap(NDView<T, 2> noise_map) {
+        noiseMap = noise_map;
+        use_noise_map = true;
+    }
+    void set_peripheralThresholdFactor(int factor) { peripheralThresholdFactor_ = factor; }
+    void find_clusters(NDView<T, 2> img);
+    void find_clusters_X(NDView<T, 2> img);
+    void rec_FillHit(int clusterIndex, int i, int j);
+    void single_pass(NDView<T, 2> img);
+    void first_pass();
+    void second_pass();
+    void store_clusters();
+
+    std::vector<Hit> steal_hits() {
+        std::vector<Hit> tmp;
+        std::swap(tmp, hits);
+        return tmp;
+    };
+    void clear_hits() { hits.clear(); };
+
+    void print_connections() {
+        fmt::print("Connections:\n");
+        for (auto it = child.begin(); it != child.end(); ++it) {
+            fmt::print("{} -> {}\n", it->first, it->second);
+        }
+    }
+    size_t total_clusters() const {
+        // TODO! fix for stealing
+        return hits.size();
+    }
+
+  private:
+    void add_link(int from, int to) {
+        // we want to add key from -> value to
+        // fmt::print("add_link({},{})\n", from, to);
+        auto it = child.find(from);
+        if (it == child.end()) {
+            child[from] = to;
+        } else {
+            // found need to disambiguate
+            if (it->second == to)
+                return;
+            else {
+                if (it->second > to) {
+                    // child[from] = to;
+                    auto old = it->second;
+                    it->second = to;
+                    add_link(old, to);
+                } else {
+                    // found value is smaller than what we want to link
+                    add_link(to, it->second);
+                }
+            }
+        }
+    }
+};
+template <typename T> int VarClusterFinder<T>::check_neighbours(int i, int j) {
+    std::vector<int> neighbour_labels;
+
+    for (int k = 0; k < 4; ++k) {
+        const auto row = i + di[k];
+        const auto col = j + dj[k];
+        if (row >= 0 && col >= 0 && row < shape_[0] && col < shape_[1]) {
+            auto tmp = labeled_.value(i + di[k], j + dj[k]);
+            if (tmp != 0)
+                neighbour_labels.push_back(tmp);
+        }
+    }
+
+    if (neighbour_labels.size() == 0) {
+        return 0;
+    } else {
+
+        // need to sort and add to union field
+        std::sort(neighbour_labels.rbegin(), neighbour_labels.rend());
+        auto first = neighbour_labels.begin();
+        auto last = std::unique(first, neighbour_labels.end());
+        if (last - first == 1)
+            return *neighbour_labels.begin();
+
+        for (auto current = first; current != last - 1; ++current) {
+            auto next = current + 1;
+            add_link(*current, *next);
+        }
+        return neighbour_labels.back(); // already sorted
+    }
+}
+
+template <typename T> void VarClusterFinder<T>::find_clusters(NDView<T, 2> img) {
+    original_ = img;
+    labeled_ = 0;
+    peripheral_labeled_ = 0;
+    current_label = 0;
+    child.clear();
+    first_pass();
+    // print_connections();
+    second_pass();
+    store_clusters();
+}
+
+template <typename T> void VarClusterFinder<T>::find_clusters_X(NDView<T, 2> img) {
+    original_ = img;
+    int clusterIndex = 0;
+    for (int i = 0; i < shape_[0]; ++i) {
+        for (int j = 0; j < shape_[1]; ++j) {
+            if (use_noise_map)
+                threshold_ = 5 * noiseMap(i, j);
+            if (original_(i, j) > threshold_) {
+                // printf("========== Cluster index: %d\n", clusterIndex);
+                rec_FillHit(clusterIndex, i, j);
+                clusterIndex++;
+            }
+        }
+    }
+    for (const auto &h : h_size)
+        hits.push_back(h.second);
+    h_size.clear();
+}
+
+template <typename T> void VarClusterFinder<T>::rec_FillHit(int clusterIndex, int i, int j) {
+    // printf("original_(%d, %d)=%f\n", i, j, original_(i,j));
+    // printf("h_size[%d].size=%d\n", clusterIndex, h_size[clusterIndex].size);
+    if (h_size[clusterIndex].size < MAX_CLUSTER_SIZE) {
+        h_size[clusterIndex].rows[h_size[clusterIndex].size] = i;
+        h_size[clusterIndex].cols[h_size[clusterIndex].size] = j;
+        h_size[clusterIndex].enes[h_size[clusterIndex].size] = original_(i, j);
+    }
+    h_size[clusterIndex].size += 1;
+    h_size[clusterIndex].energy += original_(i, j);
+    if (h_size[clusterIndex].max < original_(i, j)) {
+        h_size[clusterIndex].row = i;
+        h_size[clusterIndex].col = j;
+        h_size[clusterIndex].max = original_(i, j);
+    }
+    original_(i, j) = 0;
+
+    for (int k = 0; k < 8; ++k) { // 8 for 8-neighbour
+        const auto row = i + di_[k];
+        const auto col = j + dj_[k];
+        if (row >= 0 && col >= 0 && row < shape_[0] && col < shape_[1]) {
+            if (use_noise_map)
+                threshold_ = peripheralThresholdFactor_ * noiseMap(row, col);
+            if (original_(row, col) > threshold_) {
+                rec_FillHit(clusterIndex, row, col);
+            } else {
+                // if (h_size[clusterIndex].size < MAX_CLUSTER_SIZE){
+                //     h_size[clusterIndex].size += 1;
+                //     h_size[clusterIndex].rows[h_size[clusterIndex].size] = row;
+                //     h_size[clusterIndex].cols[h_size[clusterIndex].size] = col;
+                //     h_size[clusterIndex].enes[h_size[clusterIndex].size] = original_(row, col);
+                // }// ? weather to include peripheral pixels
+                original_(row, col) = 0; // remove peripheral pixels, to avoid potential influence for pedestal updating
+            }
+        }
+    }
+}
+
+template <typename T> void VarClusterFinder<T>::single_pass(NDView<T, 2> img) {
+    original_ = img;
+    labeled_ = 0;
+    current_label = 0;
+    child.clear();
+    first_pass();
+    // print_connections();
+    // second_pass();
+    // store_clusters();
+}
+
+template <typename T> void VarClusterFinder<T>::first_pass() {
+
+    for (int i = 0; i < original_.size(); ++i) {
+        if (use_noise_map)
+            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?
+            if (binary_(i, j)) {
+                auto tmp = check_neighbours(i, j);
+                if (tmp != 0) {
+                    labeled_(i, j) = tmp;
+                } else {
+                    labeled_(i, j) = ++current_label;
+                }
+            }
+        }
+    }
+}
+
+template <typename T> void VarClusterFinder<T>::second_pass() {
+
+    for (int64_t i = 0; i != labeled_.size(); ++i) {
+        auto current_label = labeled_(i);
+        if (current_label != 0) {
+            auto it = child.find(current_label);
+            while (it != child.end()) {
+                current_label = it->second;
+                it = child.find(current_label);
+                // do this once before doing the second pass?
+                // all values point to the final one...
+            }
+            labeled_(i) = current_label;
+        }
+    }
+}
+
+template <typename T> void VarClusterFinder<T>::store_clusters() {
+
+    // Accumulate hit information in a map
+    // Do we always have monotonic increasing
+    // labels? Then vector?
+    // here the translation is label -> Hit
+    std::unordered_map<int, Hit> h_size;
+    for (int i = 0; i < shape_[0]; ++i) {
+        for (int j = 0; j < shape_[1]; ++j) {
+            if (labeled_(i, j) != 0 || false
+                // (i-1 >= 0 and labeled_(i-1, j) != 0) or // another circle of peripheral pixels
+                // (j-1 >= 0 and labeled_(i, j-1) != 0) or
+                // (i+1 < shape_[0] and labeled_(i+1, j) != 0) or
+                // (j+1 < shape_[1] and labeled_(i, j+1) != 0)
+            ) {
+                Hit &record = h_size[labeled_(i, j)];
+                if (record.size < MAX_CLUSTER_SIZE) {
+                    record.rows[record.size] = i;
+                    record.cols[record.size] = j;
+                    record.enes[record.size] = original_(i, j);
+                } else {
+                    continue;
+                }
+                record.size += 1;
+                record.energy += original_(i, j);
+
+                if (record.max < original_(i, j)) {
+                    record.row = i;
+                    record.col = j;
+                    record.max = original_(i, j);
+                }
+            }
+        }
+    }
+
+    for (const auto &h : h_size)
+        hits.push_back(h.second);
+}
+
+} // namespace aare

include/aare/defs.hpp

@@ -0,0 +1,156 @@
+#pragma once
+
+#include "aare/Dtype.hpp"
+// #include "aare/utils/logger.hpp"
+
+#include <array>
+#include <stdexcept>
+
+#include <cassert>
+#include <cstdint>
+#include <cstring>
+#include <string>
+#include <string_view>
+#include <variant>
+#include <vector>
+
+/**
+ * @brief LOCATION macro to get the current location in the code
+ */
+#define LOCATION std::string(__FILE__) + std::string(":") + std::to_string(__LINE__) + ":" + std::string(__func__) + ":"
+
+
+namespace aare {
+
+class Cluster {
+  public:
+    int cluster_sizeX;
+    int cluster_sizeY;
+    int16_t x;
+    int16_t y;
+    Dtype dt;
+
+  private:
+    std::byte *m_data;
+
+  public:
+    Cluster(int cluster_sizeX_, int cluster_sizeY_, Dtype dt_ = Dtype(typeid(int32_t)))
+        : cluster_sizeX(cluster_sizeX_), cluster_sizeY(cluster_sizeY_), dt(dt_) {
+        m_data = new std::byte[cluster_sizeX * cluster_sizeY * dt.bytes()]{};
+    }
+    Cluster() : Cluster(3, 3) {}
+    Cluster(const Cluster &other) : Cluster(other.cluster_sizeX, other.cluster_sizeY, other.dt) {
+        if (this == &other)
+            return;
+        x = other.x;
+        y = other.y;
+        memcpy(m_data, other.m_data, other.bytes());
+    }
+    Cluster &operator=(const Cluster &other) {
+        if (this == &other)
+            return *this;
+        this->~Cluster();
+        new (this) Cluster(other);
+        return *this;
+    }
+    Cluster(Cluster &&other) noexcept
+        : cluster_sizeX(other.cluster_sizeX), cluster_sizeY(other.cluster_sizeY), x(other.x), y(other.y), dt(other.dt),
+          m_data(other.m_data) {
+        other.m_data = nullptr;
+        other.dt = Dtype(Dtype::TypeIndex::ERROR);
+    }
+    ~Cluster() { delete[] m_data; }
+    template <typename T> T get(int idx) {
+        (sizeof(T) == dt.bytes()) ? 0 : throw std::invalid_argument("[ERROR] Type size mismatch");
+        return *reinterpret_cast<T *>(m_data + idx * dt.bytes());
+    }
+    template <typename T> auto set(int idx, T val) {
+        (sizeof(T) == dt.bytes()) ? 0 : throw std::invalid_argument("[ERROR] Type size mismatch");
+        return memcpy(m_data + idx * dt.bytes(), &val, (size_t)dt.bytes());
+    }
+    // auto x() const { return x; }
+    // auto y() const { return y; }
+    // auto x(int16_t x_) { return x = x_; }
+    // auto y(int16_t y_) { return y = y_; }
+
+    template <typename T> std::string to_string() const {
+        (sizeof(T) == dt.bytes()) ? 0 : throw std::invalid_argument("[ERROR] Type size mismatch");
+        std::string s = "x: " + std::to_string(x) + " y: " + std::to_string(y) + "\nm_data: [";
+        for (int i = 0; i < cluster_sizeX * cluster_sizeY; i++) {
+            s += std::to_string(*reinterpret_cast<T *>(m_data + i * dt.bytes())) + " ";
+        }
+        s += "]";
+        return s;
+    }
+    /**
+     * @brief size of the cluster in bytes when saved to a file
+     */
+    size_t size() const { return cluster_sizeX * cluster_sizeY ; }
+    size_t bytes() const { return cluster_sizeX * cluster_sizeY * dt.bytes(); }
+    auto begin() const { return m_data; }
+    auto end() const { return m_data + cluster_sizeX * cluster_sizeY * dt.bytes(); }
+    std::byte *data() { return m_data; }
+};
+
+/**
+ * @brief header contained in parts of frames
+ */
+struct sls_detector_header {
+    uint64_t frameNumber;
+    uint32_t expLength;
+    uint32_t packetNumber;
+    uint64_t bunchId;
+    uint64_t timestamp;
+    uint16_t modId;
+    uint16_t row;
+    uint16_t column;
+    uint16_t reserved;
+    uint32_t debug;
+    uint16_t roundRNumber;
+    uint8_t detType;
+    uint8_t version;
+    std::array<uint8_t, 64> packetMask;
+    std::string to_string() {
+        std::string packetMaskStr = "[";
+        for (auto &i : packetMask) {
+            packetMaskStr += std::to_string(i) + ", ";
+        }
+        packetMaskStr += "]";
+
+        return "frameNumber: " + std::to_string(frameNumber) + "\n" + "expLength: " + std::to_string(expLength) + "\n" +
+               "packetNumber: " + std::to_string(packetNumber) + "\n" + "bunchId: " + std::to_string(bunchId) + "\n" +
+               "timestamp: " + std::to_string(timestamp) + "\n" + "modId: " + std::to_string(modId) + "\n" +
+               "row: " + std::to_string(row) + "\n" + "column: " + std::to_string(column) + "\n" +
+               "reserved: " + std::to_string(reserved) + "\n" + "debug: " + std::to_string(debug) + "\n" +
+               "roundRNumber: " + std::to_string(roundRNumber) + "\n" + "detType: " + std::to_string(detType) + "\n" +
+               "version: " + std::to_string(version) + "\n" + "packetMask: " + packetMaskStr + "\n";
+    }
+};
+
+template <typename T> struct t_xy {
+    T row;
+    T col;
+    bool operator==(const t_xy &other) const { return row == other.row && col == other.col; }
+    bool operator!=(const t_xy &other) const { return !(*this == other); }
+    std::string to_string() const { return "{ x: " + std::to_string(row) + " y: " + std::to_string(col) + " }"; }
+};
+using xy = t_xy<uint32_t>;
+
+using dynamic_shape = std::vector<int64_t>;
+
+enum class DetectorType { Jungfrau, Eiger, Mythen3, Moench, ChipTestBoard, Unknown };
+
+enum class TimingMode { Auto, Trigger };
+
+template <class T> T StringTo(const std::string &arg) { return T(arg); }
+
+template <class T> std::string toString(T arg) { return T(arg); }
+
+template <> DetectorType StringTo(const std::string & /*name*/);
+template <> std::string toString(DetectorType arg);
+
+template <> TimingMode StringTo(const std::string & /*mode*/);
+
+using DataTypeVariants = std::variant<uint16_t, uint32_t>;
+
+} // namespace aare

include/aare/json.hpp

@@ -0,0 +1,68 @@
+#pragma once
+#include <array>
+#include <map>
+#include <string>
+#include <vector>
+
+// helper functions to write json
+// append to string for better performance (not tested)
+
+namespace aare {
+
+/**
+ * @brief write a digit to a string
+ * takes key and value and outputs->"key": value,
+ * @tparam T type of value (int, uint32_t, ...)
+ * @param s string to append to
+ * @param key key to write
+ * @param value value to write
+ * @return void
+ * @note
+ * - can't use concepts here because we are using c++17
+ */
+template <typename T> inline void write_digit(std::string &s, const std::string &key, const T &value) {
+    s += "\"";
+    s += key;
+    s += "\": ";
+    s += std::to_string(value);
+    s += ", ";
+}
+inline void write_str(std::string &s, const std::string &key, const std::string &value) {
+    s += "\"";
+    s += key;
+    s += "\": \"";
+    s += value;
+    s += "\", ";
+}
+inline void write_map(std::string &s, const std::string &key, const std::map<std::string, std::string> &value) {
+    s += "\"";
+    s += key;
+    s += "\": {";
+    for (const auto &kv : value) {
+        write_str(s, kv.first, kv.second);
+    }
+    // remove last comma or trailing spaces
+    for (size_t i = s.size() - 1; i > 0; i--) {
+        if ((s[i] == ',') || (s[i] == ' ')) {
+            s.pop_back();
+        } else
+            break;
+    }
+    s += "}, ";
+}
+
+template <typename T, int N> void write_array(std::string &s, const std::string &key, const std::array<T, N> &value) {
+    s += "\"";
+    s += key;
+    s += "\": [";
+
+    for (size_t i = 0; i < N - 1; i++) {
+        s += std::to_string(value[i]);
+        s += ", ";
+    }
+    s += std::to_string(value[N - 1]);
+
+    s += "], ";
+}
+
+} // namespace aare

python/CMakeLists.txt

@@ -0,0 +1,34 @@
+if(AARE_FETCH_PYBIND11)
+    FetchContent_Declare(
+        pybind11
+        GIT_REPOSITORY https://github.com/pybind/pybind11
+        GIT_TAG        v2.11.0
+    )
+    FetchContent_MakeAvailable(pybind11)
+else()
+    find_package(pybind11 2.11 REQUIRED)
+endif()
+
+
+pybind11_add_module(
+    _aare 
+    src/module.cpp
+    )
+set_target_properties(_aare PROPERTIES
+    LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}
+)
+target_link_libraries(_aare PRIVATE aare_core aare_compiler_flags)
+
+set( PYTHON_FILES
+    aare/__init__.py
+)
+
+foreach(FILE ${PYTHON_FILES})
+    configure_file( ${FILE}  
+    ${CMAKE_BINARY_DIR}/${FILE}  )
+endforeach(FILE ${PYTHON_FILES})
+
+
+configure_file( examples/play.py  
+    ${CMAKE_BINARY_DIR}/play.py 
+)

python/aare/__init__.py

@@ -0,0 +1,3 @@
+
+
+from _aare import File

python/examples/play.py

@@ -0,0 +1,14 @@
+import matplotlib.pyplot as plt
+import numpy as np
+plt.ion()
+
+import aare
+from pathlib import Path
+
+p = Path('/Users/erik/data/aare_test_data/jungfrau/jungfrau_single_master_0.json')
+
+f = aare.File(p)
+frame = f.read_frame()
+
+fig, ax = plt.subplots()
+im = ax.imshow(frame, cmap='viridis')

python/src/file.hpp

@@ -0,0 +1,126 @@
+#include "aare/Frame.hpp"
+#include "aare/File.hpp"
+#include "aare/defs.hpp"
+// #include "aare/fClusterFileV2.hpp"
+
+#include <cstdint>
+#include <filesystem>
+#include <pybind11/numpy.h>
+#include <pybind11/iostream.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl/filesystem.h>
+#include <string>
+
+namespace py = pybind11;
+using namespace::aare;
+
+void define_file_io_bindings(py::module &m) {
+    py::class_<xy>(m, "xy")
+        .def(py::init<>())
+        .def(py::init<uint32_t, uint32_t>())
+        .def_readwrite("row", &xy::row)
+        .def_readwrite("col", &xy::col)
+        .def("__eq__", &xy::operator==)
+        .def("__ne__", &xy::operator!=)
+        .def("__repr__",
+             [](const xy &a) { return "<xy: row=" + std::to_string(a.row) + ", col=" + std::to_string(a.col) + ">"; });
+
+
+    py::class_<File>(m, "File")
+        .def(py::init([](const std::filesystem::path &fname) { return File(fname, "r", {}); }))
+        .def(
+            py::init([](const std::filesystem::path &fname, const std::string &mode) { return File(fname, mode, {}); }))
+        .def(py::init<const std::filesystem::path &, const std::string &, const FileConfig &>())
+        // .def("read", py::overload_cast<>(&File::read))
+        // .def("read", py::overload_cast<size_t>(&File::read))
+        .def("iread", py::overload_cast<size_t>(&File::iread),py::call_guard<py::gil_scoped_release>())
+        .def("frame_number", &File::frame_number)
+        .def_property_readonly("bytes_per_frame", &File::bytes_per_frame)
+        .def_property_readonly("pixels_per_frame", &File::pixels_per_frame)
+        .def("seek", &File::seek)
+        .def("tell", &File::tell)
+        .def_property_readonly("total_frames", &File::total_frames)
+        .def_property_readonly("rows", &File::rows)
+        .def_property_readonly("cols", &File::cols)
+        .def_property_readonly("bitdepth", &File::bitdepth)
+        .def_property_readonly("detector_type", &File::detector_type)
+        .def_property_readonly("geometry", &File::geometry,
+                               py::call_guard<py::scoped_ostream_redirect, py::scoped_estream_redirect>())
+        // .def("set_total_frames", &File::set_total_frames)
+        .def("read_frame", [](File &self) {
+            const uint8_t item_size = self.bytes_per_pixel();
+            py::array image;
+            std::vector<ssize_t> shape;
+            shape.reserve(2);
+            shape.push_back(self.rows());
+            shape.push_back(self.cols());
+            if (item_size == 1) {
+                image = py::array_t<uint8_t>(shape);
+            } else if (item_size == 2) {
+                image = py::array_t<uint16_t>(shape);
+            } else if (item_size == 4) {
+                image = py::array_t<uint32_t>(shape);
+            }
+            self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()));
+            return image;
+        });
+
+    py::class_<FileConfig>(m, "FileConfig")
+        .def(py::init<>())
+        .def_readwrite("rows", &FileConfig::rows)
+        .def_readwrite("cols", &FileConfig::cols)
+        .def_readwrite("version", &FileConfig::version)
+        .def_readwrite("geometry", &FileConfig::geometry)
+        .def_readwrite("detector_type", &FileConfig::detector_type)
+        .def_readwrite("max_frames_per_file", &FileConfig::max_frames_per_file)
+        .def_readwrite("total_frames", &FileConfig::total_frames)
+        .def_readwrite("dtype", &FileConfig::dtype)
+        .def("__eq__", &FileConfig::operator==)
+        .def("__ne__", &FileConfig::operator!=)
+        .def("__repr__", [](const FileConfig &a) { return "<FileConfig: " + a.to_string() + ">"; });
+
+    // py::class_<ClusterHeader>(m, "ClusterHeader")
+    //     .def(py::init<>())
+    //     .def_readwrite("frame_number", &ClusterHeader::frame_number)
+    //     .def_readwrite("n_clusters", &ClusterHeader::n_clusters)
+    //     .def("__repr__", [](const ClusterHeader &a) { return "<ClusterHeader: " + a.to_string() + ">"; });
+
+    // py::class_<ClusterV2_>(m, "ClusterV2_")
+    //     .def(py::init<>())
+    //     .def_readwrite("x", &ClusterV2_::x)
+    //     .def_readwrite("y", &ClusterV2_::y)
+    //     .def_readwrite("data", &ClusterV2_::data)
+    //     .def("__repr__", [](const ClusterV2_ &a) { return "<ClusterV2_: " + a.to_string(false) + ">"; });
+
+    // py::class_<ClusterV2>(m, "ClusterV2")
+    //     .def(py::init<>())
+    //     .def_readwrite("cluster", &ClusterV2::cluster)
+    //     .def_readwrite("frame_number", &ClusterV2::frame_number)
+    //     .def("__repr__", [](const ClusterV2 &a) { return "<ClusterV2: " + a.to_string() + ">"; });
+
+    // py::class_<ClusterFileV2>(m, "ClusterFileV2")
+    //     .def(py::init<const std::filesystem::path &, const std::string &>())
+    //     .def("read", py::overload_cast<>(&ClusterFileV2::read))
+    //     .def("read", py::overload_cast<int>(&ClusterFileV2::read))
+    //     .def("frame_number", &ClusterFileV2::frame_number)
+    //     .def("write", py::overload_cast<std::vector<ClusterV2> const &>(&ClusterFileV2::write))
+        
+    //     .def("close", &ClusterFileV2::close);
+
+    // m.def("to_clustV2", [](std::vector<Cluster> &clusters, const int frame_number) {
+    //     std::vector<ClusterV2> clusters_;
+    //     for (auto &c : clusters) {
+    //         ClusterV2 cluster;
+    //         cluster.cluster.x = c.x;
+    //         cluster.cluster.y = c.y;
+    //         int i=0;
+    //         for(auto &d : cluster.cluster.data) {
+    //             d=c.get<double>(i++);
+    //         }
+    //         cluster.frame_number = frame_number;
+    //         clusters_.push_back(cluster);
+    //     }
+    //     return clusters_;
+    // });
+}

python/src/module.cpp

@@ -0,0 +1,14 @@
+
+#include "file.hpp"
+
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+namespace py = pybind11;
+
+PYBIND11_MODULE(_aare, m) {
+
+    
+
+    define_file_io_bindings(m);
+}

python/src/np_helper.hpp

@@ -0,0 +1,115 @@
+#pragma once
+
+#include <iostream>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+#include "aare/Frame.hpp"
+#include "aare/NDArray.hpp"
+
+namespace py = pybind11;
+
+// Pass image data back to python as a numpy array
+// template <typename T, ssize_t Ndim>
+// py::array return_image_data(pl::ImageData<T, Ndim> *image) {
+
+//     py::capsule free_when_done(image, [](void *f) {
+//         pl::ImageData<T, Ndim> *foo =
+//             reinterpret_cast<pl::ImageData<T, Ndim> *>(f);
+//         delete foo;
+//     });
+
+//     return py::array_t<T>(
+//         image->shape(),        // shape
+//         image->byte_strides(), // C-style contiguous strides for double
+//         image->data(),         // the data pointer
+//         free_when_done);       // numpy array references this parent
+// }
+
+// template <typename T> py::array return_vector(std::vector<T> *vec) {
+//     py::capsule free_when_done(vec, [](void *f) {
+//         std::vector<T> *foo = reinterpret_cast<std::vector<T> *>(f);
+//         delete foo;
+//     });
+//     return py::array_t<T>({vec->size()}, // shape
+//                           {sizeof(T)}, // C-style contiguous strides for double
+//                           vec->data(), // the data pointer
+//                           free_when_done); // numpy array references this parent
+// }
+
+// template <typename Reader> py::array do_read(Reader &r, size_t n_frames) {
+//     py::array image;
+//     if (n_frames == 0)
+//         n_frames = r.total_frames();
+
+//     std::array<ssize_t, 3> shape{static_cast<ssize_t>(n_frames), r.rows(),
+//                                  r.cols()};
+//     const uint8_t item_size = r.bytes_per_pixel();
+//     if (item_size == 1) {
+//         image = py::array_t<uint8_t, py::array::c_style | py::array::forcecast>(
+//             shape);
+//     } else if (item_size == 2) {
+//         image =
+//             py::array_t<uint16_t, py::array::c_style | py::array::forcecast>(
+//                 shape);
+//     } else if (item_size == 4) {
+//         image =
+//             py::array_t<uint32_t, py::array::c_style | py::array::forcecast>(
+//                 shape);
+//     }
+//     r.read_into(reinterpret_cast<std::byte *>(image.mutable_data()), n_frames);
+//     return image;
+// }
+
+py::array return_frame(pl::Frame *ptr) {
+    py::capsule free_when_done(ptr, [](void *f) {
+        pl::Frame *foo = reinterpret_cast<pl::Frame *>(f);
+        delete foo;
+    });
+
+    const uint8_t item_size = ptr->bytes_per_pixel();
+    std::vector<ssize_t> shape;
+    for (auto val : ptr->shape())
+        if (val > 1)
+            shape.push_back(val);
+
+    std::vector<ssize_t> strides;
+    if (shape.size() == 1)
+        strides.push_back(item_size);
+    else if (shape.size() == 2) {
+        strides.push_back(item_size * shape[1]);
+        strides.push_back(item_size);
+    }
+
+    if (item_size == 1)
+        return py::array_t<uint8_t>(
+            shape, strides,
+            reinterpret_cast<uint8_t *>(ptr->data()), free_when_done);
+    else if (item_size == 2)
+        return py::array_t<uint16_t>(shape, strides,
+                                     reinterpret_cast<uint16_t *>(ptr->data()),
+                                     free_when_done);
+    else if (item_size == 4)
+        return py::array_t<uint32_t>(shape, strides,
+                                     reinterpret_cast<uint32_t *>(ptr->data()),
+                                     free_when_done);
+    return {};
+}
+
+// todo rewrite generic
+template <class T, int Flags> auto get_shape_3d(py::array_t<T, Flags> arr) {
+    return pl::Shape<3>{arr.shape(0), arr.shape(1), arr.shape(2)};
+}
+
+template <class T, int Flags> auto make_span_3d(py::array_t<T, Flags> arr) {
+    return pl::DataSpan<T, 3>(arr.mutable_data(), get_shape_3d<T, Flags>(arr));
+}
+
+template <class T, int Flags> auto get_shape_2d(py::array_t<T, Flags> arr) {
+    return pl::Shape<2>{arr.shape(0), arr.shape(1)};
+}
+
+template <class T, int Flags> auto make_span_2d(py::array_t<T, Flags> arr) {
+    return pl::DataSpan<T, 2>(arr.mutable_data(), get_shape_2d<T, Flags>(arr));
+}

src/ClusterFinder.test.cpp

@@ -0,0 +1,59 @@
+#include "aare/ClusterFinder.hpp"
+#include "aare/Pedestal.hpp"
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <chrono>
+#include <random>
+
+using namespace aare;
+
+class ClusterFinderUnitTest : public ClusterFinder {
+  public:
+    ClusterFinderUnitTest(int cluster_sizeX, int cluster_sizeY, double nSigma = 5.0, double threshold = 0.0)
+        : ClusterFinder(cluster_sizeX, cluster_sizeY, nSigma, threshold) {}
+    double get_c2() { return c2; }
+    double get_c3() { return c3; }
+    auto get_threshold() { return m_threshold; }
+    auto get_nSigma() { return m_nSigma; }
+    auto get_cluster_sizeX() { return m_cluster_sizeX; }
+    auto get_cluster_sizeY() { return m_cluster_sizeY; }
+};
+
+TEST_CASE("test ClusterFinder constructor") {
+    ClusterFinderUnitTest cf(55, 100);
+    REQUIRE(cf.get_cluster_sizeX() == 55);
+    REQUIRE(cf.get_cluster_sizeY() == 100);
+    REQUIRE(cf.get_threshold() == 0.0);
+    REQUIRE(cf.get_nSigma() == 5.0);
+    double c2 = sqrt((100 + 1) / 2 * (55 + 1) / 2);
+    double c3 = sqrt(55 * 100);
+    // REQUIRE(compare_floats<double>(cf.get_c2(), c2));
+    // REQUIRE(compare_floats<double>(cf.get_c3(), c3));
+    REQUIRE_THAT(cf.get_c2(), Catch::Matchers::WithinRel(c2, 1e-9));
+    REQUIRE_THAT(cf.get_c3(), Catch::Matchers::WithinRel(c3, 1e-9));
+}
+
+TEST_CASE("test cluster finder") {
+    aare::Pedestal pedestal(10, 10, 5);
+    NDArray<double, 2> frame({10, 10});
+    frame = 0;
+    ClusterFinder clusterFinder(3, 3, 1, 1); // 3x3 cluster, 1 nSigma, 1 threshold
+
+    auto clusters = clusterFinder.find_clusters_without_threshold(frame.span(), pedestal);
+
+    REQUIRE(clusters.size() == 0);
+
+    frame(5, 5) = 10;
+    clusters = clusterFinder.find_clusters_without_threshold(frame.span(), pedestal);
+    REQUIRE(clusters.size() == 1);
+    REQUIRE(clusters[0].x == 5);
+    REQUIRE(clusters[0].y == 5);
+    for (int i = 0; i < 3; i++) {
+        for (int j = 0; j < 3; j++) {
+            if (i == 1 && j == 1)
+                REQUIRE(clusters[0].get<double>(i * 3 + j) == 10);
+            else
+                REQUIRE(clusters[0].get<double>(i * 3 + j) == 0);
+        }
+    }
+}

src/Dtype.cpp

@@ -0,0 +1,191 @@
+
+#include "aare/Dtype.hpp"
+#include "aare/defs.hpp"
+
+#include <fmt/core.h>
+
+namespace aare {
+
+/**
+ * @brief Construct a DType object from a type_info object
+ * @param t type_info object
+ * @throw runtime_error if the type is not supported
+ * @note supported types are: int8_t, uint8_t, int16_t, uint16_t, int32_t, uint32_t, int64_t, uint64_t, float, double
+ * @note the type_info object is obtained using typeid (e.g. typeid(int))
+ */
+Dtype::Dtype(const std::type_info &t) {
+    if (t == typeid(int8_t))
+        m_type = TypeIndex::INT8;
+    else if (t == typeid(uint8_t))
+        m_type = TypeIndex::UINT8;
+    else if (t == typeid(int16_t))
+        m_type = TypeIndex::INT16;
+    else if (t == typeid(uint16_t))
+        m_type = TypeIndex::UINT16;
+    else if (t == typeid(int32_t))
+        m_type = TypeIndex::INT32;
+    else if (t == typeid(uint32_t))
+        m_type = TypeIndex::UINT32;
+    else if (t == typeid(int64_t)) // NOLINT
+        m_type = TypeIndex::INT64;
+    else if (t == typeid(uint64_t))
+        m_type = TypeIndex::UINT64;
+    else if (t == typeid(float))
+        m_type = TypeIndex::FLOAT;
+    else if (t == typeid(double))
+        m_type = TypeIndex::DOUBLE;
+    else
+        throw std::runtime_error("Could not construct data type. Type not supported.");
+}
+
+/**
+ * @brief Get the bitdepth of the data type
+ * @return bitdepth
+ */
+uint8_t Dtype::bitdepth() const {
+    switch (m_type) {
+    case TypeIndex::INT8:
+    case TypeIndex::UINT8:
+        return 8;
+    case TypeIndex::INT16:
+    case TypeIndex::UINT16:
+        return 16;
+    case TypeIndex::INT32:
+    case TypeIndex::UINT32:
+        return 32;
+    case TypeIndex::INT64:
+    case TypeIndex::UINT64:
+        return 64;
+    case TypeIndex::FLOAT:
+        return 32;
+    case TypeIndex::DOUBLE:
+        return 64;
+    case TypeIndex::NONE:
+        return 0;
+    default:
+        throw std::runtime_error(LOCATION + "Could not get bitdepth. Type not supported.");
+    }
+}
+
+/**
+ * @brief Get the number of bytes of the data type
+ */
+size_t Dtype::bytes() const { return bitdepth() / 8; }
+
+/**
+ * @brief Construct a DType object from a TypeIndex
+ * @param ti TypeIndex
+ *
+ */
+Dtype::Dtype(Dtype::TypeIndex ti) : m_type(ti) {}
+
+/**
+ * @brief Construct a DType object from a string
+ * @param sv string_view
+ * @throw runtime_error if the type is not supported
+ * @note example strings: "<i4", "u8", "f4"
+ * @note the endianess is checked and only native endianess is supported
+ */
+Dtype::Dtype(std::string_view sv) {
+
+    // Check if the file is using our native endianess
+    if (auto pos = sv.find_first_of("<>"); pos != std::string_view::npos) {
+        const auto endianess = [](const char c) {
+            if (c == '<')
+                return endian::little;
+            return endian::big;
+        }(sv[pos]);
+        if (endianess != endian::native) {
+            throw std::runtime_error("Non native endianess not supported");
+        }
+    }
+
+    // we are done with the endianess so we can remove the prefix
+    sv.remove_prefix(std::min(sv.find_first_not_of("<>"), sv.size()));
+
+    if (sv == "i1")
+        m_type = TypeIndex::INT8;
+    else if (sv == "u1")
+        m_type = TypeIndex::UINT8;
+    else if (sv == "i2")
+        m_type = TypeIndex::INT16;
+    else if (sv == "u2")
+        m_type = TypeIndex::UINT16;
+    else if (sv == "i4")
+        m_type = TypeIndex::INT32;
+    else if (sv == "u4")
+        m_type = TypeIndex::UINT32;
+    else if (sv == "i8")
+        m_type = TypeIndex::INT64;
+    else if (sv == "u8")
+        m_type = TypeIndex::UINT64;
+    else if (sv == "f4")
+        m_type = TypeIndex::FLOAT;
+    else if (sv == "f8")
+        m_type = TypeIndex::DOUBLE;
+    else
+        throw std::runtime_error("Cannot construct data type from string.");
+}
+
+Dtype Dtype::from_bitdepth(uint8_t bitdepth) {
+    switch (bitdepth) {
+    case 8:
+        return Dtype(TypeIndex::UINT8);
+    case 16:
+        return Dtype(TypeIndex::UINT16);
+    case 32:
+        return Dtype(TypeIndex::UINT32);
+    case 64:
+        return Dtype(TypeIndex::UINT64);
+    default:
+        throw std::runtime_error("Could not construct data type from bitdepth.");
+    }
+}
+/**
+ * @brief Get the string representation of the data type
+ * @return string representation
+ */
+std::string Dtype::to_string() const {
+
+    char ec{};
+    if (endian::native == endian::little)
+        ec = '<';
+    else
+        ec = '>';
+
+    switch (m_type) {
+    case TypeIndex::INT8:
+        return fmt::format("{}i1", ec);
+    case TypeIndex::UINT8:
+        return fmt::format("{}u1", ec);
+    case TypeIndex::INT16:
+        return fmt::format("{}i2", ec);
+    case TypeIndex::UINT16:
+        return fmt::format("{}u2", ec);
+    case TypeIndex::INT32:
+        return fmt::format("{}i4", ec);
+    case TypeIndex::UINT32:
+        return fmt::format("{}u4", ec);
+    case TypeIndex::INT64:
+        return fmt::format("{}i8", ec);
+    case TypeIndex::UINT64:
+        return fmt::format("{}u8", ec);
+    case TypeIndex::FLOAT:
+        return "f4";
+    case TypeIndex::DOUBLE:
+        return "f8";
+    case TypeIndex::ERROR:
+        throw std::runtime_error("Could not get string representation. Type not supported.");
+    case TypeIndex::NONE:
+        throw std::runtime_error("Could not get string representation. Type not supported.");
+    }
+    return {};
+}
+
+bool Dtype::operator==(const Dtype &other) const noexcept { return m_type == other.m_type; }
+bool Dtype::operator!=(const Dtype &other) const noexcept { return !(*this == other); }
+
+bool Dtype::operator==(const std::type_info &t) const { return Dtype(t) == *this; }
+bool Dtype::operator!=(const std::type_info &t) const { return Dtype(t) != *this; }
+
+} // namespace aare

src/Dtype.test.cpp

@@ -0,0 +1,54 @@
+
+
+#include "aare/Dtype.hpp"
+#include <catch2/catch_test_macros.hpp>
+
+using aare::Dtype;
+using aare::endian;
+
+TEST_CASE("Construct from typeid") {
+    REQUIRE(Dtype(typeid(int)) == typeid(int));
+    REQUIRE(Dtype(typeid(int)) != typeid(double));
+}
+
+TEST_CASE("Construct from string") {
+    if (endian::native == endian::little) {
+        REQUIRE(Dtype("<i1") == typeid(int8_t));
+        REQUIRE(Dtype("<u1") == typeid(uint8_t));
+        REQUIRE(Dtype("<i2") == typeid(int16_t));
+        REQUIRE(Dtype("<u2") == typeid(uint16_t));
+        REQUIRE(Dtype("<i4") == typeid(int));
+        REQUIRE(Dtype("<u4") == typeid(unsigned));
+        REQUIRE(Dtype("<i4") == typeid(int32_t));
+        // REQUIRE(Dtype("<i8") == typeid(long));
+        REQUIRE(Dtype("<i8") == typeid(int64_t));
+        REQUIRE(Dtype("<u4") == typeid(uint32_t));
+        REQUIRE(Dtype("<u8") == typeid(uint64_t));
+        REQUIRE(Dtype("f4") == typeid(float));
+        REQUIRE(Dtype("f8") == typeid(double));
+    }
+
+    if (endian::native == endian::big) {
+        REQUIRE(Dtype(">i1") == typeid(int8_t));
+        REQUIRE(Dtype(">u1") == typeid(uint8_t));
+        REQUIRE(Dtype(">i2") == typeid(int16_t));
+        REQUIRE(Dtype(">u2") == typeid(uint16_t));
+        REQUIRE(Dtype(">i4") == typeid(int));
+        REQUIRE(Dtype(">u4") == typeid(unsigned));
+        REQUIRE(Dtype(">i4") == typeid(int32_t));
+        // REQUIRE(Dtype(">i8") == typeid(long));
+        REQUIRE(Dtype(">i8") == typeid(int64_t));
+        REQUIRE(Dtype(">u4") == typeid(uint32_t));
+        REQUIRE(Dtype(">u8") == typeid(uint64_t));
+        REQUIRE(Dtype("f4") == typeid(float));
+        REQUIRE(Dtype("f8") == typeid(double));
+    }
+}
+
+TEST_CASE("Construct from string with endianess") {
+    // TODO! handle big endian system in test!
+    REQUIRE(Dtype("<i4") == typeid(int32_t));
+    REQUIRE_THROWS(Dtype(">i4") == typeid(int32_t));
+}
+
+TEST_CASE("Convert to string") { REQUIRE(Dtype(typeid(int)).to_string() == "<i4"); }

src/File.cpp

@@ -0,0 +1,71 @@
+#include "aare/File.hpp"
+#include "aare/NumpyFile.hpp"
+#include "aare/RawFile.hpp"
+
+#include <fmt/format.h>
+
+namespace aare {
+
+File::File(const std::filesystem::path &fname, const std::string &mode, const FileConfig &cfg)
+    : file_impl(nullptr), is_npy(true) {
+    if (mode != "r" && mode != "w" && mode != "a") {
+        throw std::invalid_argument("Unsupported file mode");
+    }
+
+    if ((mode == "r" || mode == "a") && !std::filesystem::exists(fname)) {
+        throw std::runtime_error(fmt::format("File does not exist: {}", fname.string()));
+    }
+
+    if (fname.extension() == ".raw" || fname.extension() == ".json") {
+        // aare::logger::debug("Loading raw file");
+        file_impl = new RawFile(fname, mode, cfg);
+        is_npy = false;
+    }
+    // check if extension is numpy
+    else if (fname.extension() == ".npy") {
+        // aare::logger::debug("Loading numpy file");
+        file_impl = new NumpyFile(fname, mode, cfg);
+    } else {
+        throw std::runtime_error("Unsupported file type");
+    }
+}
+
+void File::write(Frame &frame, sls_detector_header header) {
+    if (is_npy) {
+        // aare::logger::info("ignoring header for npy file");
+        dynamic_cast<NumpyFile *>(file_impl)->write(frame);
+    } else {
+        dynamic_cast<RawFile *>(file_impl)->write(frame, header);
+    }
+}
+Frame File::read() { return file_impl->read(); }
+size_t File::total_frames() const { return file_impl->total_frames(); }
+std::vector<Frame> File::read(size_t n_frames) { return file_impl->read(n_frames); }
+void File::read_into(std::byte *image_buf) { file_impl->read_into(image_buf); }
+void File::read_into(std::byte *image_buf, size_t n_frames) { file_impl->read_into(image_buf, n_frames); }
+size_t File::frame_number(size_t frame_index) { return file_impl->frame_number(frame_index); }
+size_t File::bytes_per_frame() { return file_impl->bytes_per_frame(); }
+size_t File::pixels_per_frame() { return file_impl->pixels_per_frame(); }
+void File::seek(size_t frame_number) { file_impl->seek(frame_number); }
+size_t File::tell() const { return file_impl->tell(); }
+size_t File::rows() const { return file_impl->rows(); }
+size_t File::cols() const { return file_impl->cols(); }
+size_t File::bitdepth() const { return file_impl->bitdepth(); }
+size_t File::bytes_per_pixel() const { return file_impl->bitdepth()/8; }
+void File::set_total_frames(size_t total_frames) { return file_impl->set_total_frames(total_frames); }
+File::~File() { delete file_impl; }
+DetectorType File::detector_type() const { return file_impl->detector_type(); }
+xy File::geometry() const {
+    if (is_npy) {
+        return {1, 1};
+    }
+    return reinterpret_cast<RawFile *>(file_impl)->geometry();
+}
+
+Frame File::iread(size_t frame_number) { return file_impl->iread(frame_number); }
+
+File::File(File &&other) noexcept : file_impl(other.file_impl), is_npy(other.is_npy) { other.file_impl = nullptr; }
+
+// write move assignment operator
+
+} // namespace aare

src/Frame.cpp

@@ -0,0 +1,110 @@
+#include "aare/Frame.hpp"
+
+#include <cstddef>
+#include <cstring>
+#include <iostream>
+#include <sys/types.h>
+
+namespace aare {
+
+/**
+ * @brief Construct a new Frame
+ * @param bytes pointer to the data to be copied into the frame
+ * @param rows number of rows
+ * @param cols number of columns
+ * @param bitdepth bitdepth of the pixels
+ */
+Frame::Frame(const std::byte *bytes, uint32_t rows, uint32_t cols, Dtype dtype)
+    : m_rows(rows), m_cols(cols), m_dtype(dtype), m_data(new std::byte[rows * cols * m_dtype.bytes()]) {
+
+    std::memcpy(m_data, bytes, rows * cols * m_dtype.bytes());
+}
+
+/**
+ * @brief Construct a new Frame
+ * @param rows number of rows
+ * @param cols number of columns
+ * @param bitdepth bitdepth of the pixels
+ * @note the data is initialized to zero
+ */
+Frame::Frame(uint32_t rows, uint32_t cols, Dtype dtype)
+    : m_rows(rows), m_cols(cols), m_dtype(dtype), m_data(new std::byte[rows * cols * dtype.bytes()]) {
+
+    std::memset(m_data, 0, rows * cols * dtype.bytes());
+}
+
+uint32_t Frame::rows() const { return m_rows; }
+uint32_t Frame::cols() const { return m_cols; }
+size_t Frame::bitdepth() const { return m_dtype.bitdepth(); }
+Dtype Frame::dtype() const { return m_dtype; }
+uint64_t Frame::size() const { return m_rows * m_cols; }
+size_t Frame::bytes() const { return m_rows * m_cols * m_dtype.bytes(); }
+std::byte *Frame::data() const { return m_data; }
+
+/**
+ * @brief Get the pointer to the pixel at the given row and column
+ * @param row row index
+ * @param col column index
+ * @return pointer to the pixel
+ * @note the user should cast the pointer to the appropriate type
+ */
+std::byte *Frame::get(uint32_t row, uint32_t col) {
+    if ((row >= m_rows) || (col >= m_cols)) {
+        std::cerr << "Invalid row or column index" << '\n';
+        return nullptr;
+    }
+    return m_data + (row * m_cols + col) * (m_dtype.bytes());
+}
+
+// Frame &Frame::operator=(const Frame &other) {
+//     if (this == &other) {
+//         return *this;
+//     }
+//     m_rows = other.rows();
+//     m_cols = other.cols();
+//     m_dtype = other.dtype();
+//     m_data = new std::byte[m_rows * m_cols * m_dtype.bytes()];
+//     if (m_data == nullptr) {
+//         throw std::bad_alloc();
+//     }
+//     std::memcpy(m_data, other.m_data, m_rows * m_cols * m_dtype.bytes());
+//     return *this;
+// }
+Frame &Frame::operator=(Frame &&other) noexcept {
+    if (this == &other) {
+        return *this;
+    }
+    m_rows = other.rows();
+    m_cols = other.cols();
+    m_dtype = other.dtype();
+    if (m_data != nullptr) {
+        delete[] m_data;
+    }
+    m_data = other.m_data;
+    other.m_data = nullptr;
+    other.m_rows = other.m_cols = 0;
+    other.m_dtype = Dtype(Dtype::TypeIndex::ERROR);
+    return *this;
+}
+Frame::Frame(Frame &&other) noexcept
+    : m_rows(other.rows()), m_cols(other.cols()), m_dtype(other.dtype()), m_data(other.m_data) {
+
+    other.m_data = nullptr;
+    other.m_rows = other.m_cols = 0;
+    other.m_dtype = Dtype(Dtype::TypeIndex::ERROR);
+}
+// Frame::Frame(const Frame &other)
+//     : m_rows(other.rows()), m_cols(other.cols()), m_dtype(other.dtype()),
+//       m_data(new std::byte[m_rows * m_cols * m_dtype.bytes()]) {
+
+//     std::memcpy(m_data, other.m_data, m_rows * m_cols * m_dtype.bytes());
+// }
+
+Frame Frame::copy() const {
+    Frame frame(m_rows, m_cols, m_dtype);
+    std::memcpy(frame.m_data, m_data, m_rows * m_cols * m_dtype.bytes());
+    return frame;
+}
+
+Frame::~Frame() noexcept { delete[] m_data; }
+} // namespace aare

src/Frame.test.cpp

@@ -0,0 +1,152 @@
+#include "aare/Frame.hpp"
+#include "aare/Dtype.hpp"
+#include <catch2/catch_test_macros.hpp>
+
+using namespace aare;
+
+TEST_CASE("Construct a frame") {
+    size_t rows = 10;
+    size_t cols = 10;
+    size_t bitdepth = 8;
+
+    Frame frame(rows, cols, Dtype::from_bitdepth(bitdepth));
+
+    REQUIRE(frame.rows() == rows);
+    REQUIRE(frame.cols() == cols);
+    REQUIRE(frame.bitdepth() == bitdepth);
+    REQUIRE(frame.bytes() == rows * cols * bitdepth / 8);
+
+    // data should be initialized to 0
+    for (size_t i = 0; i < rows; i++) {
+        for (size_t j = 0; j < cols; j++) {
+            uint8_t *data = (uint8_t *)frame.get(i, j);
+            REQUIRE(data != nullptr);
+            REQUIRE(*data == 0);
+        }
+    }
+}
+
+TEST_CASE("Set a value in a 8 bit frame") {
+    size_t rows = 10;
+    size_t cols = 10;
+    size_t bitdepth = 8;
+
+    Frame frame(rows, cols, Dtype::from_bitdepth(bitdepth));
+
+    // set a value
+    uint8_t value = 255;
+    frame.set(5, 7, value);
+
+    // only the value we did set should be non-zero
+    for (size_t i = 0; i < rows; i++) {
+        for (size_t j = 0; j < cols; j++) {
+            uint8_t *data = (uint8_t *)frame.get(i, j);
+            REQUIRE(data != nullptr);
+            if (i == 5 && j == 7) {
+                REQUIRE(*data == value);
+            } else {
+                REQUIRE(*data == 0);
+            }
+        }
+    }
+}
+
+TEST_CASE("Set a value in a 64 bit frame") {
+    size_t rows = 10;
+    size_t cols = 10;
+    size_t bitdepth = 64;
+
+    Frame frame(rows, cols, Dtype::from_bitdepth(bitdepth));
+
+    // set a value
+    uint64_t value = 255;
+    frame.set(5, 7, value);
+
+    // only the value we did set should be non-zero
+    for (size_t i = 0; i < rows; i++) {
+        for (size_t j = 0; j < cols; j++) {
+            uint64_t *data = (uint64_t *)frame.get(i, j);
+            REQUIRE(data != nullptr);
+            if (i == 5 && j == 7) {
+                REQUIRE(*data == value);
+            } else {
+                REQUIRE(*data == 0);
+            }
+        }
+    }
+}
+
+TEST_CASE("Move construct a frame") {
+    size_t rows = 10;
+    size_t cols = 10;
+    size_t bitdepth = 8;
+
+    Frame frame(rows, cols, Dtype::from_bitdepth(bitdepth));
+    std::byte *data = frame.data();
+
+    Frame frame2(std::move(frame));
+
+    // state of the moved from object
+    REQUIRE(frame.rows() == 0);
+    REQUIRE(frame.cols() == 0);
+    REQUIRE(frame.dtype() == Dtype(Dtype::TypeIndex::ERROR));
+    REQUIRE(frame.data() == nullptr);
+
+    // state of the moved to object
+    REQUIRE(frame2.rows() == rows);
+    REQUIRE(frame2.cols() == cols);
+    REQUIRE(frame2.bitdepth() == bitdepth);
+    REQUIRE(frame2.bytes() == rows * cols * bitdepth / 8);
+    REQUIRE(frame2.data() == data);
+}
+
+TEST_CASE("Move assign a frame") {
+    size_t rows = 10;
+    size_t cols = 10;
+    size_t bitdepth = 8;
+
+    Frame frame(rows, cols, Dtype::from_bitdepth(bitdepth));
+    std::byte *data = frame.data();
+
+    Frame frame2(5, 5, Dtype::from_bitdepth(16));
+
+    frame2 = std::move(frame);
+
+    // state of the moved from object
+    REQUIRE(frame.rows() == 0);
+    REQUIRE(frame.cols() == 0);
+    REQUIRE(frame.dtype() == Dtype(Dtype::TypeIndex::ERROR));
+    REQUIRE(frame.data() == nullptr);
+
+    // state of the moved to object
+    REQUIRE(frame2.rows() == rows);
+    REQUIRE(frame2.cols() == cols);
+    REQUIRE(frame2.bitdepth() == bitdepth);
+    REQUIRE(frame2.bytes() == rows * cols * bitdepth / 8);
+    REQUIRE(frame2.data() == data);
+}
+
+TEST_CASE("test explicit copy constructor") {
+    size_t rows = 10;
+    size_t cols = 10;
+    size_t bitdepth = 8;
+
+    Frame frame(rows, cols, Dtype::from_bitdepth(bitdepth));
+    std::byte *data = frame.data();
+
+    Frame frame2 = frame.copy();
+
+    // state of the original object
+    REQUIRE(frame.rows() == rows);
+    REQUIRE(frame.cols() == cols);
+    REQUIRE(frame.bitdepth() == bitdepth);
+    REQUIRE(frame.bytes() == rows * cols * bitdepth / 8);
+    REQUIRE(frame.data() == data);
+
+    // state of the copied object
+    REQUIRE(frame2.rows() == rows);
+    REQUIRE(frame2.cols() == cols);
+    REQUIRE(frame2.bitdepth() == bitdepth);
+    REQUIRE(frame2.bytes() == rows * cols * bitdepth / 8);
+    REQUIRE(frame2.data() != data);
+}

src/NDArray.test.cpp

@@ -0,0 +1,377 @@
+#include "aare/NDArray.hpp"
+#include <array>
+#include <catch2/catch_test_macros.hpp>
+
+using aare::NDArray;
+using aare::NDView;
+using aare::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 (uint32_t i = 0; i < image.size(); ++i) {
+        REQUIRE(image(i) == view(i));
+    }
+
+    // Changing the image doesn't change the view
+    image = 43;
+    for (uint32_t i = 0; i < image.size(); ++i) {
+        REQUIRE(image(i) != view(i));
+    }
+}
+
+TEST_CASE("1D image") {
+    std::array<int64_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<int64_t, 2> shape{{3, 7}};
+    NDArray<long> img(shape, 5);
+    for (uint32_t i = 0; i != img.size(); ++i) {
+        REQUIRE(img(i) == 5);
+    }
+
+    for (uint32_t 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 (uint32_t i = 0; i != img.size(); ++i) {
+        REQUIRE(img(i) == 5.0f);
+    }
+
+    // Double check general properties
+    REQUIRE(img.size() == 3 * 4 * 2);
+
+    for (uint32_t 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<int64_t, 3> shape{3, 4, 2};
+    NDArray<double, 3> a{shape};
+    NDArray<double, 3> b{shape};
+    for (uint32_t 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") {
+    int64_t w = 15;
+    int64_t h = 75;
+    std::array<int64_t, 2> shape{w, h};
+    NDArray<double> a{shape};
+    REQUIRE(a.size() == static_cast<uint64_t>(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 (uint32_t 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<int64_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<int64_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 (uint32_t i = 0; i < C.size(); ++i) {
+            REQUIRE(C(i) == a_val + b_val);
+        }
+
+        // Value of A is not changed
+        for (uint32_t i = 0; i < A.size(); ++i) {
+            REQUIRE(A(i) == a_val);
+        }
+
+        // Value of B is not changed
+        for (uint32_t 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 (uint32_t i = 0; i < C.size(); ++i) {
+            REQUIRE(C(i) == a_val - b_val);
+        }
+
+        // Value of A is not changed
+        for (uint32_t i = 0; i < A.size(); ++i) {
+            REQUIRE(A(i) == a_val);
+        }
+
+        // Value of B is not changed
+        for (uint32_t 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 (uint32_t i = 0; i < C.size(); ++i) {
+            REQUIRE(C(i) == a_val * b_val);
+        }
+
+        // Value of A is not changed
+        for (uint32_t i = 0; i < A.size(); ++i) {
+            REQUIRE(A(i) == a_val);
+        }
+
+        // Value of B is not changed
+        for (uint32_t 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 (uint32_t i = 0; i < C.size(); ++i) {
+            REQUIRE(C(i) == a_val / b_val);
+        }
+
+        // Value of A is not changed
+        for (uint32_t i = 0; i < A.size(); ++i) {
+            REQUIRE(A(i) == a_val);
+        }
+
+        // Value of B is not changed
+        for (uint32_t 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 (uint32_t i = 0; i < C.size(); ++i) {
+            REQUIRE(C(i) == a_val - v);
+        }
+
+        // Value of A is not changed
+        for (uint32_t 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 (uint32_t i = 0; i < C.size(); ++i) {
+            REQUIRE(C(i) == a_val + v);
+        }
+
+        // Value of A is not changed
+        for (uint32_t 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 (uint32_t i = 0; i < C.size(); ++i) {
+            REQUIRE(C(i) == a_val / v);
+        }
+
+        // Value of A is not changed
+        for (uint32_t 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 (uint32_t i = 0; i < C.size(); ++i) {
+            REQUIRE(C(i) == a_val / v);
+        }
+
+        // Value of A is not changed
+        for (uint32_t i = 0; i < A.size(); ++i) {
+            REQUIRE(A(i) == a_val);
+        }
+    }
+}

src/NDView.test.cpp

@@ -0,0 +1,193 @@
+#include "aare/NDView.hpp"
+#include <catch2/catch_test_macros.hpp>
+
+#include <iostream>
+#include <vector>
+
+using aare::NDView;
+using aare::Shape;
+
+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<int64_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<int64_t>(vec0.size())});
+    NDView<int, 1> data1(vec1.data(), Shape<1>{static_cast<int64_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);
+}
+
+TEST_CASE("compare two views") {
+    std::vector<int> vec1;
+    for (int i = 0; i != 12; ++i) {
+        vec1.push_back(i);
+    }
+    NDView<int, 2> view1(vec1.data(), Shape<2>{3, 4});
+
+    std::vector<int> vec2;
+    for (int i = 0; i != 12; ++i) {
+        vec2.push_back(i);
+    }
+    NDView<int, 2> view2(vec2.data(), Shape<2>{3, 4});
+
+    REQUIRE((view1 == view2));
+}

src/NumpyFile.cpp

@@ -0,0 +1,200 @@
+
+#include "aare/NumpyFile.hpp"
+#include "aare/NumpyHelpers.hpp"
+
+namespace aare {
+
+
+
+NumpyFile::NumpyFile(const std::filesystem::path &fname, const std::string &mode, FileConfig cfg) {
+    // TODO! add opts to constructor
+    m_fname = fname;
+    m_mode = mode;
+    if (mode == "r") {
+        fp = fopen(m_fname.string().c_str(), "rb");
+        if (!fp) {
+            throw std::runtime_error(fmt::format("Could not open: {} for reading", m_fname.string()));
+        }
+        load_metadata();
+    } else if (mode == "w") {
+        m_bitdepth = cfg.dtype.bitdepth();
+        m_rows = cfg.rows;
+        m_cols = cfg.cols;
+        m_header = {cfg.dtype, false, {cfg.rows, cfg.cols}};
+        m_header.shape = {0, cfg.rows, cfg.cols};
+        fp = fopen(m_fname.string().c_str(), "wb");
+        if (!fp) {
+            throw std::runtime_error(fmt::format("Could not open: {} for reading", m_fname.string()));
+        }
+        initial_header_len = aare::NumpyHelpers::write_header(std::filesystem::path(m_fname.c_str()), m_header);
+    }
+    m_pixels_per_frame = std::accumulate(m_header.shape.begin() + 1, m_header.shape.end(), 1, std::multiplies<>());
+
+    m_bytes_per_frame = m_header.dtype.bitdepth() / 8 * m_pixels_per_frame;
+}
+void NumpyFile::write(Frame &frame) { write_impl(frame.data(), frame.bytes()); }
+void NumpyFile::write_impl(void *data, uint64_t size) {
+
+    if (fp == nullptr) {
+        throw std::runtime_error("File not open");
+    }
+    if (!(m_mode == "w" || m_mode == "a")) {
+        throw std::invalid_argument("File not open for writing");
+    }
+    if (fseek(fp, 0, SEEK_END))
+        throw std::runtime_error("Could not seek to end of file");
+    size_t const rc = fwrite(data, size, 1, fp);
+    if (rc != 1) {
+        throw std::runtime_error("Error writing frame to file");
+    }
+
+    m_header.shape[0]++;
+}
+
+Frame NumpyFile::get_frame(size_t frame_number) {
+    Frame frame(m_header.shape[1], m_header.shape[2], m_header.dtype);
+    get_frame_into(frame_number, frame.data());
+    return frame;
+}
+void NumpyFile::get_frame_into(size_t frame_number, std::byte *image_buf) {
+    if (fp == nullptr) {
+        throw std::runtime_error("File not open");
+    }
+    if (frame_number > m_header.shape[0]) {
+        throw std::invalid_argument("Frame number out of range");
+    }
+    if (fseek(fp, header_size + frame_number * m_bytes_per_frame, SEEK_SET)) // NOLINT
+        throw std::runtime_error("Could not seek to frame");
+
+    size_t const rc = fread(image_buf, m_bytes_per_frame, 1, fp);
+    if (rc != 1) {
+        throw std::runtime_error("Error reading frame from file");
+    }
+}
+
+size_t NumpyFile::pixels_per_frame() { return m_pixels_per_frame; };
+size_t NumpyFile::bytes_per_frame() { return m_bytes_per_frame; };
+
+std::vector<Frame> NumpyFile::read(size_t n_frames) {
+    // TODO: implement this in a more efficient way
+    std::vector<Frame> frames;
+    for (size_t i = 0; i < n_frames; i++) {
+        frames.push_back(get_frame(current_frame));
+        current_frame++;
+    }
+    return frames;
+}
+void NumpyFile::read_into(std::byte *image_buf, size_t n_frames) {
+    // TODO: implement this in a more efficient way
+    for (size_t i = 0; i < n_frames; i++) {
+        get_frame_into(current_frame++, image_buf);
+        image_buf += m_bytes_per_frame;
+    }
+}
+
+NumpyFile::~NumpyFile() noexcept {
+    if (m_mode == "w" || m_mode == "a") {
+        // determine number of frames
+        if (fseek(fp, 0, SEEK_END)) {
+            std::cout << "Could not seek to end of file" << std::endl;
+        }
+        size_t const file_size = ftell(fp);
+        size_t const data_size = file_size - initial_header_len;
+        size_t const n_frames = data_size / m_bytes_per_frame;
+        // update number of frames in header (first element of shape)
+        m_header.shape[0] = n_frames;
+        if (fseek(fp, 0, SEEK_SET)) {
+            std::cout << "Could not seek to beginning of file" << std::endl;
+        }
+        // create string stream to contain header
+        std::stringstream ss;
+        aare::NumpyHelpers::write_header(ss, m_header);
+        std::string const header_str = ss.str();
+        // write header
+        size_t const rc = fwrite(header_str.c_str(), header_str.size(), 1, fp);
+        if (rc != 1) {
+            std::cout << "Error writing header to numpy file in destructor" << std::endl;
+        }
+    }
+
+    if (fp != nullptr) {
+        if (fclose(fp)) {
+            std::cout << "Error closing file" << std::endl;
+        }
+    }
+}
+
+void NumpyFile::load_metadata() {
+
+    // read magic number
+    std::array<char, 6> tmp{};
+    size_t rc = fread(tmp.data(), tmp.size(), 1, fp);
+    if (rc != 1) {
+        throw std::runtime_error("Error reading magic number");
+    }
+    if (tmp != aare::NumpyHelpers::magic_str) {
+        for (auto item : tmp)
+            fmt::print("{}, ", static_cast<int>(item));
+        fmt::print("\n");
+        throw std::runtime_error("Not a numpy file");
+    }
+
+    // read version
+    rc = fread(reinterpret_cast<char *>(&major_ver_), sizeof(major_ver_), 1, fp);
+    rc += fread(reinterpret_cast<char *>(&minor_ver_), sizeof(minor_ver_), 1, fp);
+    if (rc != 2) {
+        throw std::runtime_error("Error reading numpy version");
+    }
+
+    if (major_ver_ == 1) {
+        header_len_size = 2;
+    } else if (major_ver_ == 2) {
+        header_len_size = 4;
+    } else {
+        throw std::runtime_error("Unsupported numpy version");
+    }
+
+    // read header length
+    rc = fread(reinterpret_cast<char *>(&header_len), header_len_size, 1, fp);
+    if (rc != 1) {
+        throw std::runtime_error("Error reading header length");
+    }
+    header_size = aare::NumpyHelpers::magic_string_length + 2 + header_len_size + header_len;
+    if (header_size % 16 != 0) {
+        fmt::print("Warning: header length is not a multiple of 16\n");
+    }
+
+    // read header
+    std::string header(header_len, '\0');
+    rc = fread(header.data(), header_len, 1, fp);
+    if (rc != 1) {
+        throw std::runtime_error("Error reading header");
+    }
+
+    // parse header
+    std::vector<std::string> const keys{"descr", "fortran_order", "shape"};
+    auto dict_map = aare::NumpyHelpers::parse_dict(header, keys);
+    if (dict_map.empty())
+        throw std::runtime_error("invalid dictionary in header");
+
+    std::string const descr_s = dict_map["descr"];
+    std::string const fortran_s = dict_map["fortran_order"];
+    std::string const shape_s = dict_map["shape"];
+
+    std::string const descr = aare::NumpyHelpers::parse_str(descr_s);
+    aare::Dtype const dtype = aare::NumpyHelpers::parse_descr(descr);
+
+    // convert literal Python bool to C++ bool
+    bool const fortran_order = aare::NumpyHelpers::parse_bool(fortran_s);
+
+    // parse the shape tuple
+    auto shape_v = aare::NumpyHelpers::parse_tuple(shape_s);
+    std::vector<size_t> shape;
+    for (const auto &item : shape_v) {
+        auto dim = static_cast<size_t>(std::stoul(item));
+        shape.push_back(dim);
+    }
+    m_header = {dtype, fortran_order, shape};
+}
+
+} // namespace aare

src/NumpyFile.test.cpp

@@ -0,0 +1,50 @@
+#include "aare/NumpyFile.hpp"
+#include "aare/NDArray.hpp"
+
+#include <catch2/catch_test_macros.hpp>
+#include "test_config.hpp"
+
+using aare::Dtype;
+using aare::NumpyFile;
+TEST_CASE("Read a 1D numpy file with int32 data type") {
+
+    auto fpath = test_data_path() / "numpy" / "test_1d_int32.npy";
+    REQUIRE(std::filesystem::exists(fpath));
+
+    NumpyFile f(fpath);
+
+    // we know the file contains 10 elements of np.int32 containing values 0-9
+    REQUIRE(f.dtype() == Dtype::INT32);
+    REQUIRE(f.shape() == std::vector<size_t>{10});
+
+    // use the load function to read the full file into a NDArray
+    auto data = f.load<int32_t, 1>();
+    for (int32_t i = 0; i < 10; i++) {
+        REQUIRE(data(i) == i);
+    }
+}
+
+TEST_CASE("Read a 3D numpy file with np.double data type") {
+
+    auto fpath = test_data_path() / "numpy" / "test_3d_double.npy";
+    REQUIRE(std::filesystem::exists(fpath));
+
+    NumpyFile f(fpath);
+
+    // we know the file contains 10 elements of np.int32 containing values 0-9
+    REQUIRE(f.dtype() == Dtype::DOUBLE);
+    REQUIRE(f.shape() == std::vector<size_t>{3, 2, 5});
+
+    // use the load function to read the full file into a NDArray
+    // numpy code to generate the array
+    //  arr2[0,0,0] = 1.0
+    //  arr2[0,0,1] = 2.0
+    //  arr2[0,1,0] = 72.0
+    //  arr2[2,0,4] = 63.0
+
+    auto data = f.load<double, 3>();
+    REQUIRE(data(0, 0, 0) == 1.0);
+    REQUIRE(data(0, 0, 1) == 2.0);
+    REQUIRE(data(0, 1, 0) == 72.0);
+    REQUIRE(data(2, 0, 4) == 63.0);
+}

src/NumpyHelpers.cpp

@@ -0,0 +1,269 @@
+/*
+   28-03-2024 modified by: Bechir Braham <bechir.braham@psi.ch>
+
+   Copyright 2017-2023 Leon Merten Lohse
+
+   Permission is hereby granted, free of charge, to any person obtaining a copy
+   of this software and associated documentation files (the "Software"), to deal
+   in the Software without restriction, including without limitation the rights
+   to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+   copies of the Software, and to permit persons to whom the Software is
+   furnished to do so, subject to the following conditions:
+
+   The above copyright notice and this permission notice shall be included in
+   all copies or substantial portions of the Software.
+
+   THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+   OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+   SOFTWARE.
+*/
+
+#include "aare/NumpyHelpers.hpp"
+#include <iterator>
+
+namespace aare {
+
+std::string NumpyHeader::to_string() const {
+    std::stringstream sstm;
+    sstm << "dtype: " << dtype.to_string() << ", fortran_order: " << fortran_order << ' ';
+    sstm << "shape: (";
+    for (auto item : shape)
+        sstm << item << ',';
+    sstm << ')';
+    return sstm.str();
+}
+
+
+namespace NumpyHelpers {
+
+std::unordered_map<std::string, std::string> parse_dict(std::string in, const std::vector<std::string> &keys) {
+    std::unordered_map<std::string, std::string> map;
+    if (keys.empty())
+        return map;
+
+    in = trim(in);
+
+    // unwrap dictionary
+    if ((in.front() == '{') && (in.back() == '}'))
+        in = in.substr(1, in.length() - 2);
+    else
+        throw std::runtime_error("Not a Python dictionary.");
+
+    std::vector<std::pair<size_t, std::string>> positions;
+
+    for (auto const &key : keys) {
+        size_t const pos = in.find("'" + key + "'");
+
+        if (pos == std::string::npos)
+            throw std::runtime_error("Missing '" + key + "' key.");
+
+        std::pair<size_t, std::string> const position_pair{pos, key};
+        positions.push_back(position_pair);
+    }
+
+    // sort by position in dict
+    std::sort(positions.begin(), positions.end());
+
+    for (size_t i = 0; i < positions.size(); ++i) {
+        std::string raw_value;
+        size_t const begin{positions[i].first};
+        size_t end{std::string::npos};
+
+        std::string const key = positions[i].second;
+
+        if (i + 1 < positions.size())
+            end = positions[i + 1].first;
+
+        raw_value = in.substr(begin, end - begin);
+
+        raw_value = trim(raw_value);
+
+        if (raw_value.back() == ',')
+            raw_value.pop_back();
+
+        map[key] = get_value_from_map(raw_value);
+    }
+
+    return map;
+}
+
+aare::Dtype parse_descr(std::string typestring) {
+
+    if (typestring.length() < 3) {
+        throw std::runtime_error("invalid typestring (length)");
+    }
+
+    constexpr char little_endian_char = '<';
+    constexpr char big_endian_char = '>';
+    constexpr char no_endian_char = '|';
+    constexpr std::array<char, 3> endian_chars = {little_endian_char, big_endian_char, no_endian_char};
+    constexpr std::array<char, 4> numtype_chars = {'f', 'i', 'u', 'c'};
+
+    const char byteorder_c = typestring[0];
+    const char kind_c = typestring[1];
+    std::string const itemsize_s = typestring.substr(2);
+
+    if (!in_array(byteorder_c, endian_chars)) {
+        throw std::runtime_error("invalid typestring (byteorder)");
+    }
+
+    if (!in_array(kind_c, numtype_chars)) {
+        throw std::runtime_error("invalid typestring (kind)");
+    }
+
+    if (!is_digits(itemsize_s)) {
+        throw std::runtime_error("invalid typestring (itemsize)");
+    }
+
+    return aare::Dtype(typestring);
+}
+
+bool parse_bool(const std::string &in) {
+    if (in == "True")
+        return true;
+    if (in == "False")
+        return false;
+    throw std::runtime_error("Invalid python boolean.");
+}
+
+std::string get_value_from_map(const std::string &mapstr) {
+    size_t const sep_pos = mapstr.find_first_of(':');
+    if (sep_pos == std::string::npos)
+        return "";
+
+    std::string const tmp = mapstr.substr(sep_pos + 1);
+    return trim(tmp);
+}
+
+bool is_digits(const std::string &str) { return std::all_of(str.begin(), str.end(), ::isdigit); }
+
+std::vector<std::string> parse_tuple(std::string in) {
+    std::vector<std::string> v;
+    const char separator = ',';
+
+    in = trim(in);
+
+    if ((in.front() == '(') && (in.back() == ')'))
+        in = in.substr(1, in.length() - 2);
+    else
+        throw std::runtime_error("Invalid Python tuple.");
+
+    std::istringstream iss(in);
+
+    for (std::string token; std::getline(iss, token, separator);) {
+        v.push_back(token);
+    }
+
+    return v;
+}
+
+std::string trim(const std::string &str) {
+    const std::string whitespace = " \t\n";
+    auto begin = str.find_first_not_of(whitespace);
+
+    if (begin == std::string::npos)
+        return "";
+
+    auto end = str.find_last_not_of(whitespace);
+    return str.substr(begin, end - begin + 1);
+}
+
+std::string parse_str(const std::string &in) {
+    if ((in.front() == '\'') && (in.back() == '\''))
+        return in.substr(1, in.length() - 2);
+
+    throw std::runtime_error("Invalid python string.");
+}
+
+void write_magic(std::ostream &ostream, int version_major, int version_minor) {
+    ostream.write(magic_str.data(), magic_string_length);
+    ostream.put(static_cast<char>(version_major));
+    ostream.put(static_cast<char>(version_minor));
+}
+template <typename T> inline std::string write_tuple(const std::vector<T> &v) {
+
+    if (v.empty())
+        return "()";
+    std::ostringstream ss;
+    ss.imbue(std::locale("C"));
+
+    if (v.size() == 1) {
+        ss << "(" << v.front() << ",)";
+    } else {
+        const std::string delimiter = ", ";
+        // v.size() > 1
+        ss << "(";
+        // for (size_t i = 0; i < v.size() - 1; ++i) {
+        //     ss << v[i] << delimiter;
+        // }
+        // ss << v.back();
+        std::copy(v.begin(), v.end() - 1, std::ostream_iterator<T>(ss, ", "));
+        ss << v.back();
+        ss << ")";
+    }
+
+    return ss.str();
+}
+
+inline std::string write_boolean(bool b) {
+    if (b)
+        return "True";
+    return "False";
+}
+
+inline std::string write_header_dict(const std::string &descr, bool fortran_order, const std::vector<size_t> &shape) {
+    std::string const s_fortran_order = write_boolean(fortran_order);
+    std::string const shape_s = write_tuple(shape);
+
+    return "{'descr': '" + descr + "', 'fortran_order': " + s_fortran_order + ", 'shape': " + shape_s + ", }";
+}
+
+size_t write_header(const std::filesystem::path &fname, const NumpyHeader &header) {
+    std::ofstream out(fname, std::ios::binary | std::ios::out);
+    return write_header(out, header);
+}
+
+size_t write_header(std::ostream &out, const NumpyHeader &header) {
+    std::string const header_dict = write_header_dict(header.dtype.to_string(), header.fortran_order, header.shape);
+
+    size_t length = magic_string_length + 2 + 2 + header_dict.length() + 1;
+
+    int version_major = 1;
+    int version_minor = 0;
+    if (length >= static_cast<size_t>(255) * 255) {
+        length = magic_string_length + 2 + 4 + header_dict.length() + 1;
+        version_major = 2;
+        version_minor = 0;
+    }
+    size_t const padding_len = 16 - length % 16;
+    std::string const padding(padding_len, ' ');
+
+    // write magic
+    write_magic(out, version_major, version_minor);
+
+    // write header length
+    if (version_major == 1 && version_minor == 0) {
+        auto header_len = static_cast<uint16_t>(header_dict.length() + padding.length() + 1);
+
+        std::array<uint8_t, 2> header_len_le16{static_cast<uint8_t>((header_len >> 0) & 0xff),
+                                               static_cast<uint8_t>((header_len >> 8) & 0xff)};
+        out.write(reinterpret_cast<char *>(header_len_le16.data()), 2);
+    } else {
+        auto header_len = static_cast<uint32_t>(header_dict.length() + padding.length() + 1);
+
+        std::array<uint8_t, 4> header_len_le32{
+            static_cast<uint8_t>((header_len >> 0) & 0xff), static_cast<uint8_t>((header_len >> 8) & 0xff),
+            static_cast<uint8_t>((header_len >> 16) & 0xff), static_cast<uint8_t>((header_len >> 24) & 0xff)};
+        out.write(reinterpret_cast<char *>(header_len_le32.data()), 4);
+    }
+
+    out << header_dict << padding << '\n';
+    return length;
+}
+
+} // namespace NumpyHelpers
+} // namespace aare

src/NumpyHelpers.test.cpp

@@ -0,0 +1,62 @@
+#include "aare/NumpyHelpers.hpp" //Is this really a public header?
+#include <catch2/catch_test_macros.hpp>
+
+using namespace aare::NumpyHelpers;
+
+TEST_CASE("is_digits with a few standard cases") {
+    REQUIRE(is_digits(""));
+    REQUIRE(is_digits("123"));
+    REQUIRE(is_digits("0"));
+    REQUIRE_FALSE(is_digits("hej123"));
+    REQUIRE_FALSE(is_digits("a"));
+    REQUIRE_FALSE(is_digits(" "));
+    REQUIRE_FALSE(is_digits("abcdef"));
+}
+
+TEST_CASE("Check for quotes and return stripped string") {
+    REQUIRE(parse_str("'hej'") == "hej");
+    REQUIRE(parse_str("'hej hej'") == "hej hej");
+    REQUIRE(parse_str("''") == "");
+}
+
+TEST_CASE("parsing a string without quotes throws") { REQUIRE_THROWS(parse_str("hej")); }
+
+TEST_CASE("trim whitespace") {
+    REQUIRE(trim(" hej ") == "hej");
+    REQUIRE(trim("hej") == "hej");
+    REQUIRE(trim(" hej") == "hej");
+    REQUIRE(trim("hej ") == "hej");
+    REQUIRE(trim(" ") == "");
+    REQUIRE(trim(" \thej hej ") == "hej hej");
+}
+
+TEST_CASE("parse data type descriptions") {
+    REQUIRE(parse_descr("<i1") == aare::Dtype::INT8);
+    REQUIRE(parse_descr("<i2") == aare::Dtype::INT16);
+    REQUIRE(parse_descr("<i4") == aare::Dtype::INT32);
+    REQUIRE(parse_descr("<i8") == aare::Dtype::INT64);
+
+    REQUIRE(parse_descr("<u1") == aare::Dtype::UINT8);
+    REQUIRE(parse_descr("<u2") == aare::Dtype::UINT16);
+    REQUIRE(parse_descr("<u4") == aare::Dtype::UINT32);
+    REQUIRE(parse_descr("<u8") == aare::Dtype::UINT64);
+
+    REQUIRE(parse_descr("<f4") == aare::Dtype::FLOAT);
+    REQUIRE(parse_descr("<f8") == aare::Dtype::DOUBLE);
+}
+
+TEST_CASE("is element in array") {
+    REQUIRE(in_array(1, std::array<int, 3>{1, 2, 3}));
+    REQUIRE_FALSE(in_array(4, std::array<int, 3>{1, 2, 3}));
+    REQUIRE(in_array(1, std::array<int, 1>{1}));
+    REQUIRE_FALSE(in_array(1, std::array<int, 0>{}));
+}
+
+TEST_CASE("Parse numpy dict") {
+    std::string in = "{'descr': '<f4', 'fortran_order': False, 'shape': (3, 4)}";
+    std::vector<std::string> keys{"descr", "fortran_order", "shape"};
+    auto map = parse_dict(in, keys);
+    REQUIRE(map["descr"] == "'<f4'");
+    REQUIRE(map["fortran_order"] == "False");
+    REQUIRE(map["shape"] == "(3, 4)");
+}

src/Pedestal.test.cpp

@@ -0,0 +1,103 @@
+#include "aare/Pedestal.hpp"
+
+
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+#include <catch2/catch_test_macros.hpp>
+#include <chrono>
+#include <random>
+
+using namespace aare;
+TEST_CASE("test pedestal constructor") {
+    aare::Pedestal pedestal(10, 10, 5);
+    REQUIRE(pedestal.rows() == 10);
+    REQUIRE(pedestal.cols() == 10);
+    REQUIRE(pedestal.n_samples() == 5);
+    for (int i = 0; i < 10; i++) {
+        for (int j = 0; j < 10; j++) {
+            REQUIRE(pedestal.get_sum()(i, j) == 0);
+            REQUIRE(pedestal.get_sum2()(i, j) == 0);
+            REQUIRE(pedestal.cur_samples()(i, j) == 0);
+        }
+    }
+}
+
+TEST_CASE("test pedestal push") {
+    aare::Pedestal pedestal(10, 10, 5);
+    aare::Frame frame(10, 10, Dtype::UINT16);
+    for (int i = 0; i < 10; i++) {
+        for (int j = 0; j < 10; j++) {
+            frame.set<uint16_t>(i, j, i + j);
+        }
+    }
+
+    // test single push
+    pedestal.push<uint16_t>(frame);
+    for (int i = 0; i < 10; i++) {
+        for (int j = 0; j < 10; j++) {
+            REQUIRE(pedestal.get_sum()(i, j) == i + j);
+            REQUIRE(pedestal.get_sum2()(i, j) == (i + j) * (i + j));
+            REQUIRE(pedestal.cur_samples()(i, j) == 1);
+        }
+    }
+
+    // test clear
+    pedestal.clear();
+    for (int i = 0; i < 10; i++) {
+        for (int j = 0; j < 10; j++) {
+            REQUIRE(pedestal.get_sum()(i, j) == 0);
+            REQUIRE(pedestal.get_sum2()(i, j) == 0);
+            REQUIRE(pedestal.cur_samples()(i, j) == 0);
+        }
+    }
+
+    // test number of samples after multiple push
+    for (uint32_t k = 0; k < 50; k++) {
+        pedestal.push<uint16_t>(frame);
+        for (uint32_t i = 0; i < 10; i++) {
+            for (uint32_t j = 0; j < 10; j++) {
+                if (k < 5) {
+                    REQUIRE(pedestal.cur_samples()(i, j) == k + 1);
+                    REQUIRE(pedestal.get_sum()(i, j) == (k + 1) * (i + j));
+                    REQUIRE(pedestal.get_sum2()(i, j) == (k + 1) * (i + j) * (i + j));
+                } else {
+                    REQUIRE(pedestal.cur_samples()(i, j) == 5);
+                    REQUIRE(pedestal.get_sum()(i, j) == 5 * (i + j));
+                    REQUIRE(pedestal.get_sum2()(i, j) == 5 * (i + j) * (i + j));
+                }
+                REQUIRE(pedestal.mean(i, j) == (i + j));
+                REQUIRE(pedestal.variance(i, j) == 0);
+                REQUIRE(pedestal.standard_deviation(i, j) == 0);
+            }
+        }
+    }
+}
+
+TEST_CASE("test pedestal with normal distribution") {
+    const double MEAN = 5.0, STD = 2.0, VAR = STD * STD, TOLERANCE = 0.1;
+
+    unsigned seed = std::chrono::system_clock::now().time_since_epoch().count();
+    std::default_random_engine generator(seed);
+    std::normal_distribution<double> distribution(MEAN, STD);
+
+    aare::Pedestal pedestal(3, 5, 10000);
+    for (int i = 0; i < 10000; i++) {
+        aare::Frame frame(3, 5, Dtype::DOUBLE);
+        for (int j = 0; j < 3; j++) {
+            for (int k = 0; k < 5; k++) {
+                frame.set<double>(j, k, distribution(generator));
+            }
+        }
+        pedestal.push<double>(frame);
+    }
+    auto mean = pedestal.mean();
+    auto variance = pedestal.variance();
+    auto standard_deviation = pedestal.standard_deviation();
+
+    for (int i = 0; i < 3; i++) {
+        for (int j = 0; j < 5; j++) {
+            REQUIRE_THAT(mean(i, j), Catch::Matchers::WithinAbs(MEAN, MEAN * TOLERANCE));
+            REQUIRE_THAT(variance(i, j), Catch::Matchers::WithinAbs(VAR, VAR * TOLERANCE));
+            REQUIRE_THAT(standard_deviation(i, j), Catch::Matchers::WithinAbs(STD, STD * TOLERANCE));
+        }
+    }
+}

src/RawFile.cpp

@@ -0,0 +1,428 @@
+#include "aare/RawFile.hpp"
+#include "aare/defs.hpp"
+#include "aare/json.hpp"
+
+#include <fmt/format.h>
+#include <nlohmann/json.hpp>
+
+using json = nlohmann::json;
+
+namespace aare {
+
+RawFile::RawFile(const std::filesystem::path &fname, const std::string &mode, const FileConfig &config) {
+    m_mode = mode;
+    m_fname = fname;
+    if (mode == "r" || mode == "r+") {
+        if (config != FileConfig()) {
+            // aare::logger::warn(
+            //     "In read mode it is not necessary to provide a config, the provided config will be ignored");
+        }
+        parse_fname();
+        parse_metadata();
+        find_number_of_subfiles();
+        find_geometry();
+        open_subfiles();
+
+    } else if (mode == "w" || mode == "w+") {
+
+        if (std::filesystem::exists(fname)) {
+            // handle mode w as w+ (no overrwriting)
+            throw std::runtime_error(LOCATION + "File already exists");
+        }
+
+        parse_config(config);
+        parse_fname();
+        write_master_file();
+        n_subfiles = 1;
+        n_subfile_parts = 1;
+        subfile_cols = m_cols;
+        subfile_rows = m_rows;
+        open_subfiles();
+
+    } else {
+        throw std::runtime_error(LOCATION + "Unsupported mode");
+    }
+}
+
+void RawFile::parse_config(const FileConfig &config) {
+    m_bitdepth = config.dtype.bitdepth();
+    m_total_frames = config.total_frames;
+    m_rows = config.rows;
+    m_cols = config.cols;
+    m_type = config.detector_type;
+    max_frames_per_file = config.max_frames_per_file;
+    m_geometry = config.geometry;
+    version = config.version;
+    subfile_rows = config.geometry.row;
+    subfile_cols = config.geometry.col;
+
+    if (m_geometry != aare::xy{1, 1}) {
+        throw std::runtime_error(LOCATION + "Only geometry {1,1} files are supported for writing");
+    }
+}
+void RawFile::write_master_file() {
+    if (m_ext != ".json") {
+        throw std::runtime_error(LOCATION + "only json master files are supported for writing");
+    }
+    std::ofstream ofs(master_fname(), std::ios::binary);
+    std::string ss;
+    ss.reserve(1024);
+    ss += "{\n\t";
+    aare::write_str(ss, "Version", version);
+    ss += "\n\t";
+    aare::write_digit(ss, "Total Frames", m_total_frames);
+    ss += "\n\t";
+    aare::write_str(ss, "Detector Type", toString(m_type));
+    ss += "\n\t";
+    aare::write_str(ss, "Geometry", m_geometry.to_string());
+    ss += "\n\t";
+
+    uint64_t img_size = (m_cols * m_rows) / (static_cast<size_t>(m_geometry.col * m_geometry.row));
+    img_size *= m_bitdepth;
+    aare::write_digit(ss, "Image Size in bytes", img_size);
+    ss += "\n\t";
+    aare::write_digit(ss, "Max Frames Per File", max_frames_per_file);
+    ss += "\n\t";
+    aare::write_digit(ss, "Dynamic Range", m_bitdepth);
+    ss += "\n\t";
+    const aare::xy pixels = {static_cast<uint32_t>(m_rows / m_geometry.row),
+                             static_cast<uint32_t>(m_cols / m_geometry.col)};
+    aare::write_str(ss, "Pixels", pixels.to_string());
+    ss += "\n\t";
+    aare::write_digit(ss, "Number of rows", m_rows);
+    ss += "\n\t";
+    const std::string tmp = "{\n"
+                            "        \"Frame Number\": \"8 bytes\",\n"
+                            "        \"Exposure Length\": \"4 bytes\",\n"
+                            "        \"Packet Number\": \"4 bytes\",\n"
+                            "        \"Bunch Id\": \"8 bytes\",\n"
+                            "        \"Timestamp\": \"8 bytes\",\n"
+                            "        \"Module Id\": \"2 bytes\",\n"
+                            "        \"Row\": \"2 bytes\",\n"
+                            "        \"Column\": \"2 bytes\",\n"
+                            "        \"Reserved\": \"2 bytes\",\n"
+                            "        \"Debug\": \"4 bytes\",\n"
+                            "        \"RoundRNumber\": \"2 bytes\",\n"
+                            "        \"DetType\": \"1 byte\",\n"
+                            "        \"Version\": \"1 byte\",\n"
+                            "        \"Packet Mask\": \"64 bytes\"\n"
+                            "    }";
+
+    ss += "\"Frame Header Format\":" + tmp + "\n";
+    ss += "}";
+    ofs << ss;
+    ofs.close();
+}
+
+void RawFile::open_subfiles() {
+    if (m_mode == "r")
+        for (size_t i = 0; i != n_subfiles; ++i) {
+            auto v = std::vector<SubFile *>(n_subfile_parts);
+            for (size_t j = 0; j != n_subfile_parts; ++j) {
+                v[j] = new SubFile(data_fname(i, j), m_type, subfile_rows, subfile_cols, m_bitdepth);
+            }
+            subfiles.push_back(v);
+        }
+    else {
+        auto v = std::vector<SubFile *>(n_subfile_parts); // only one subfile is implemented
+        v[0] = new SubFile(data_fname(0, 0), m_type, m_rows, m_cols, m_bitdepth, "w");
+        subfiles.push_back(v);
+    }
+}
+
+sls_detector_header RawFile::read_header(const std::filesystem::path &fname) {
+    sls_detector_header h{};
+    FILE *fp = fopen(fname.string().c_str(), "r");
+    if (!fp)
+        throw std::runtime_error(fmt::format("Could not open: {} for reading", fname.string()));
+
+    size_t const rc = fread(reinterpret_cast<char *>(&h), sizeof(h), 1, fp);
+    if (rc != 1)
+        throw std::runtime_error(LOCATION + "Could not read header from file");
+    if (fclose(fp)) {
+        throw std::runtime_error(LOCATION + "Could not close file");
+    }
+
+    return h;
+}
+bool RawFile::is_master_file(const std::filesystem::path &fpath) {
+    std::string const stem = fpath.stem().string();
+    return stem.find("_master_") != std::string::npos;
+}
+
+void RawFile::find_number_of_subfiles() {
+    int n_mod = 0;
+    while (std::filesystem::exists(data_fname(++n_mod, 0)))
+        ;
+    n_subfiles = n_mod;
+}
+inline std::filesystem::path RawFile::data_fname(size_t mod_id, size_t file_id) {
+    return this->m_base_path / fmt::format("{}_d{}_f{}_{}.raw", this->m_base_name, file_id, mod_id, this->m_findex);
+}
+
+inline std::filesystem::path RawFile::master_fname() {
+    return this->m_base_path / fmt::format("{}_master_{}{}", this->m_base_name, this->m_findex, this->m_ext);
+}
+
+void RawFile::find_geometry() {
+    uint16_t r{};
+    uint16_t c{};
+    for (size_t i = 0; i < n_subfile_parts; i++) {
+        for (size_t j = 0; j != n_subfiles; ++j) {
+            auto h = this->read_header(data_fname(j, i));
+            r = std::max(r, h.row);
+            c = std::max(c, h.column);
+
+            positions.push_back({h.row, h.column});
+        }
+    }
+
+    r++;
+    c++;
+
+    m_rows = (r * subfile_rows);
+    m_cols = (c * subfile_cols);
+
+    m_rows += static_cast<size_t>((r - 1) * cfg.module_gap_row);
+}
+
+void RawFile::parse_metadata() {
+    if (m_ext == ".raw") {
+        parse_raw_metadata();
+        if (m_bitdepth == 0) {
+            switch (m_type) {
+            case DetectorType::Eiger:
+                m_bitdepth = 32;
+                break;
+            default:
+                m_bitdepth = 16;
+            }
+        }
+    } else if (m_ext == ".json") {
+        parse_json_metadata();
+    } else {
+        throw std::runtime_error(LOCATION + "Unsupported file type");
+    }
+    n_subfile_parts = static_cast<size_t>(m_geometry.row) * m_geometry.col;
+}
+
+void RawFile::parse_json_metadata() {
+    std::ifstream ifs(master_fname());
+    json j;
+    ifs >> j;
+    double v = j["Version"];
+    version = fmt::format("{:.1f}", v);
+    m_type = StringTo<DetectorType>(j["Detector Type"].get<std::string>());
+    timing_mode = StringTo<TimingMode>(j["Timing Mode"].get<std::string>());
+    m_total_frames = j["Frames in File"];
+    subfile_rows = j["Pixels"]["y"];
+    subfile_cols = j["Pixels"]["x"];
+    max_frames_per_file = j["Max Frames Per File"];
+    try {
+        m_bitdepth = j.at("Dynamic Range");
+    } catch (const json::out_of_range &e) {
+        m_bitdepth = 16;
+    }
+    // only Eiger had quad
+    if (m_type == DetectorType::Eiger) {
+        quad = (j["Quad"] == 1);
+    }
+
+    m_geometry = {j["Geometry"]["y"], j["Geometry"]["x"]};
+}
+void RawFile::parse_raw_metadata() {
+    std::ifstream ifs(master_fname());
+    for (std::string line; std::getline(ifs, line);) {
+        if (line == "#Frame Header")
+            break;
+        auto pos = line.find(':');
+        auto key_pos = pos;
+        while (key_pos != std::string::npos && std::isspace(line[--key_pos]))
+            ;
+        if (key_pos != std::string::npos) {
+            auto key = line.substr(0, key_pos + 1);
+            auto value = line.substr(pos + 2);
+            // do the actual parsing
+            if (key == "Version") {
+                version = value;
+            } else if (key == "TimeStamp") {
+
+            } else if (key == "Detector Type") {
+                m_type = StringTo<DetectorType>(value);
+            } else if (key == "Timing Mode") {
+                timing_mode = StringTo<TimingMode>(value);
+            } else if (key == "Pixels") {
+                // Total number of pixels cannot be found yet looking at
+                // submodule
+                pos = value.find(',');
+                subfile_cols = std::stoi(value.substr(1, pos));
+                subfile_rows = std::stoi(value.substr(pos + 1));
+            } else if (key == "Total Frames") {
+                m_total_frames = std::stoi(value);
+            } else if (key == "Dynamic Range") {
+                m_bitdepth = std::stoi(value);
+            } else if (key == "Quad") {
+                quad = (value == "1");
+            } else if (key == "Max Frames Per File") {
+                max_frames_per_file = std::stoi(value);
+            } else if (key == "Geometry") {
+                pos = value.find(',');
+                m_geometry = {static_cast<uint32_t>(std::stoi(value.substr(1, pos))),
+                              static_cast<uint32_t>(std::stoi(value.substr(pos + 1)))};
+            }
+        }
+    }
+}
+
+void RawFile::parse_fname() {
+    bool wrong_format = false;
+    m_base_path = m_fname.parent_path().string();
+    m_base_name = m_fname.stem().string();
+    m_ext = m_fname.extension().string();
+    try {
+        auto pos = m_base_name.rfind('_');
+        m_findex = std::stoi(m_base_name.substr(pos + 1));
+    } catch (const std::invalid_argument &e) {
+        m_findex = 0;
+        wrong_format = true;
+    }
+    auto pos = m_base_name.find("_master_");
+    if (pos != std::string::npos) {
+        m_base_name.erase(pos);
+        wrong_format = true;
+    }
+    if (wrong_format && (m_mode == "w+" || m_mode == "w")) {
+        // aare::logger::warn("Master Filename", m_fname, "is not in the correct format");
+        // aare::logger::warn("using", master_fname(), "as the master file");
+    }
+}
+
+Frame RawFile::get_frame(size_t frame_index) {
+    auto f = Frame(this->m_rows, this->m_cols, Dtype::from_bitdepth(this->m_bitdepth));
+    std::byte *frame_buffer = f.data();
+    get_frame_into(frame_index, frame_buffer);
+    return f;
+}
+
+void RawFile::get_frame_into(size_t frame_index, std::byte *frame_buffer) {
+    if (frame_index > this->m_total_frames) {
+        throw std::runtime_error(LOCATION + "Frame number out of range");
+    }
+    std::vector<size_t> frame_numbers(this->n_subfile_parts);
+    std::vector<size_t> frame_indices(this->n_subfile_parts, frame_index);
+
+    if (n_subfile_parts != 1) {
+        for (size_t part_idx = 0; part_idx != this->n_subfile_parts; ++part_idx) {
+            auto subfile_id = frame_index / this->max_frames_per_file;
+            frame_numbers[part_idx] =
+                this->subfiles[subfile_id][part_idx]->frame_number(frame_index % this->max_frames_per_file);
+        }
+        // 1. if frame number vector is the same break
+        while (std::adjacent_find(frame_numbers.begin(), frame_numbers.end(), std::not_equal_to<>()) !=
+               frame_numbers.end()) {
+            // 2. find the index of the minimum frame number,
+            auto min_frame_idx =
+                std::distance(frame_numbers.begin(), std::min_element(frame_numbers.begin(), frame_numbers.end()));
+            // 3. increase its index and update its respective frame number
+            frame_indices[min_frame_idx]++;
+            // 4. if we can't increase its index => throw error
+            if (frame_indices[min_frame_idx] >= this->m_total_frames) {
+                throw std::runtime_error(LOCATION + "Frame number out of range");
+            }
+            auto subfile_id = frame_indices[min_frame_idx] / this->max_frames_per_file;
+            frame_numbers[min_frame_idx] = this->subfiles[subfile_id][min_frame_idx]->frame_number(
+                frame_indices[min_frame_idx] % this->max_frames_per_file);
+        }
+    }
+
+    if (this->m_geometry.col == 1) {
+        // get the part from each subfile and copy it to the frame
+        for (size_t part_idx = 0; part_idx != this->n_subfile_parts; ++part_idx) {
+            auto corrected_idx = frame_indices[part_idx];
+            auto subfile_id = corrected_idx / this->max_frames_per_file;
+            auto part_offset = this->subfiles[subfile_id][part_idx]->bytes_per_part();
+            this->subfiles[subfile_id][part_idx]->get_part(frame_buffer + part_idx * part_offset,
+                                                           corrected_idx % this->max_frames_per_file);
+        }
+
+    } else {
+
+        // create a buffer that will hold a the frame part
+        auto bytes_per_part = this->subfile_rows * this->subfile_cols * this->m_bitdepth / 8;
+        auto *part_buffer = new std::byte[bytes_per_part];
+
+        for (size_t part_idx = 0; part_idx != this->n_subfile_parts; ++part_idx) {
+            auto corrected_idx = frame_indices[part_idx];
+            auto subfile_id = corrected_idx / this->max_frames_per_file;
+
+            this->subfiles[subfile_id][part_idx]->get_part(part_buffer, corrected_idx % this->max_frames_per_file);
+            for (size_t cur_row = 0; cur_row < (this->subfile_rows); cur_row++) {
+                auto irow = cur_row + (part_idx / this->m_geometry.col) * this->subfile_rows;
+                auto icol = (part_idx % this->m_geometry.col) * this->subfile_cols;
+                auto dest = (irow * this->m_cols + icol);
+                dest = dest * this->m_bitdepth / 8;
+                memcpy(frame_buffer + dest, part_buffer + cur_row * this->subfile_cols * this->m_bitdepth / 8,
+                       this->subfile_cols * this->m_bitdepth / 8);
+            }
+        }
+        delete[] part_buffer;
+    }
+}
+
+void RawFile::write(Frame &frame, sls_detector_header header) {
+    if (m_mode == "r") {
+        throw std::runtime_error(LOCATION + "File is open in read mode");
+    }
+    size_t const subfile_id = this->current_frame / this->max_frames_per_file;
+    for (size_t part_idx = 0; part_idx != this->n_subfile_parts; ++part_idx) {
+
+        this->subfiles[subfile_id][part_idx]->write_part(frame.data(), header,
+                                                         this->current_frame % this->max_frames_per_file);
+    }
+    this->current_frame++;
+}
+
+std::vector<Frame> RawFile::read(size_t n_frames) {
+    // TODO: implement this in a more efficient way
+    std::vector<Frame> frames;
+    for (size_t i = 0; i < n_frames; i++) {
+        frames.push_back(this->get_frame(this->current_frame));
+        this->current_frame++;
+    }
+    return frames;
+}
+void RawFile::read_into(std::byte *image_buf, size_t n_frames) {
+    // TODO: implement this in a more efficient way
+    for (size_t i = 0; i < n_frames; i++) {
+        this->get_frame_into(this->current_frame++, image_buf);
+        image_buf += this->bytes_per_frame();
+    }
+}
+
+size_t RawFile::frame_number(size_t frame_index) {
+    if (frame_index > this->m_total_frames) {
+        throw std::runtime_error(LOCATION + "Frame number out of range");
+    }
+    size_t const subfile_id = frame_index / this->max_frames_per_file;
+    return this->subfiles[subfile_id][0]->frame_number(frame_index % this->max_frames_per_file);
+}
+
+RawFile::~RawFile() noexcept {
+
+    // update master file
+    if (m_mode == "w" || m_mode == "w+" || m_mode == "r+") {
+        try {
+            write_master_file();
+        } catch (...) {
+            // aare::logger::warn(LOCATION + "Could not update master file");
+        }
+    }
+
+    for (auto &vec : subfiles) {
+        for (auto *subfile : vec) {
+            delete subfile;
+        }
+    }
+}
+
+} // namespace aare

src/RawFile.test.cpp

@@ -0,0 +1,114 @@
+#include "aare/File.hpp"
+
+#include <catch2/catch_test_macros.hpp>
+#include <filesystem>
+
+#include "test_config.hpp"
+
+using aare::File;
+
+TEST_CASE("Read number of frames from a jungfrau raw file") {
+
+    auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json";
+    REQUIRE(std::filesystem::exists(fpath));
+
+    File f(fpath, "r");
+    REQUIRE(f.total_frames() == 10);
+}
+
+TEST_CASE("Read frame numbers from a jungfrau raw file") {
+    auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json";
+    REQUIRE(std::filesystem::exists(fpath));
+
+    File f(fpath, "r");
+
+    // we know this file has 10 frames with frame numbers 1 to 10
+    // f0 1,2,3
+    // f1 4,5,6
+    // f2 7,8,9
+    // f3 10
+    for (size_t i = 0; i < 10; i++) {
+        CHECK(f.frame_number(i) == i + 1);
+    }
+}
+
+TEST_CASE("Read data from a jungfrau 500k single port raw file") {
+    auto fpath = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json";
+    REQUIRE(std::filesystem::exists(fpath));
+
+    File f(fpath, "r");
+
+    // we know this file has 10 frames with pixel 0,0 being: 2123, 2051, 2109, 2117, 2089, 2095, 2072, 2126, 2097, 2102
+    std::vector<uint16_t> pixel_0_0 = {2123, 2051, 2109, 2117, 2089, 2095, 2072, 2126, 2097, 2102};
+    for (size_t i = 0; i < 10; i++) {
+        auto frame = f.read();
+        CHECK(frame.rows() == 512);
+        CHECK(frame.cols() == 1024);
+        CHECK(frame.view<uint16_t>()(0, 0) == pixel_0_0[i]);
+    }
+}
+
+TEST_CASE("Read frame numbers from a raw file") {
+    auto fpath = test_data_path() / "eiger" / "eiger_500k_16bit_master_0.json";
+    REQUIRE(std::filesystem::exists(fpath));
+
+    // we know this file has 3 frames with frame numbers 14, 15, 16
+    std::vector<size_t> frame_numbers = {14, 15, 16};
+
+    File f(fpath, "r");
+    for (size_t i = 0; i < 3; i++) {
+        CHECK(f.frame_number(i) == frame_numbers[i]);
+    }
+}
+
+TEST_CASE("Compare reading from a numpy file with a raw file") {
+    auto fpath_raw = test_data_path() / "jungfrau" / "jungfrau_single_master_0.json";
+    REQUIRE(std::filesystem::exists(fpath_raw));
+
+    auto fpath_npy = test_data_path() / "jungfrau" / "jungfrau_single_0.npy";
+    REQUIRE(std::filesystem::exists(fpath_npy));
+
+    File raw(fpath_raw, "r");
+    File npy(fpath_npy, "r");
+
+    CHECK(raw.total_frames() == 10);
+    CHECK(npy.total_frames() == 10);
+
+    for (size_t i = 0; i < 10; ++i) {
+        auto raw_frame = raw.read();
+        auto npy_frame = npy.read();
+        CHECK((raw_frame.view<uint16_t>() == npy_frame.view<uint16_t>()));
+    }
+}
+
+TEST_CASE("Read multipart files") {
+    auto fpath = test_data_path() / "jungfrau" / "jungfrau_double_master_0.json";
+    REQUIRE(std::filesystem::exists(fpath));
+
+    File f(fpath, "r");
+
+    // we know this file has 10 frames check read_multiport.py for the values
+    std::vector<uint16_t> pixel_0_0 = {2099, 2121, 2108, 2084, 2084, 2118, 2066, 2108, 2112, 2116};
+    std::vector<uint16_t> pixel_0_1 = {2842, 2796, 2865, 2798, 2805, 2817, 2852, 2789, 2792, 2833};
+    std::vector<uint16_t> pixel_255_1023 = {2149, 2037, 2115, 2102, 2118, 2090, 2036, 2071, 2073, 2142};
+    std::vector<uint16_t> pixel_511_1023 = {3231, 3169, 3167, 3162, 3168, 3160, 3171, 3171, 3169, 3171};
+    std::vector<uint16_t> pixel_1_0 = {2748, 2614, 2665, 2629, 2618, 2630, 2631, 2634, 2577, 2598};
+
+    for (size_t i = 0; i < 10; i++) {
+        auto frame = f.read();
+        CHECK(frame.rows() == 512);
+        CHECK(frame.cols() == 1024);
+        CHECK(frame.view<uint16_t>()(0, 0) == pixel_0_0[i]);
+        CHECK(frame.view<uint16_t>()(0, 1) == pixel_0_1[i]);
+        CHECK(frame.view<uint16_t>()(1, 0) == pixel_1_0[i]);
+        CHECK(frame.view<uint16_t>()(255, 1023) == pixel_255_1023[i]);
+        CHECK(frame.view<uint16_t>()(511, 1023) == pixel_511_1023[i]);
+    }
+}
+
+TEST_CASE("Read file with unordered frames") {
+    auto fpath = test_data_path() / "mythen" / "scan242_master_3.raw";
+    REQUIRE(std::filesystem::exists(fpath));
+    File f(fpath, "r");
+    REQUIRE_THROWS((f.read()));
+}

src/SubFile.cpp

@@ -0,0 +1,71 @@
+#include "aare/SubFile.hpp"
+
+#include <cstring> // memcpy
+#include <fmt/core.h>
+#include <iostream>
+
+
+namespace aare {
+
+SubFile::SubFile(const std::filesystem::path &fname, DetectorType detector, size_t rows, size_t cols, size_t bitdepth,
+                 const std::string &mode)
+    : m_bitdepth(bitdepth), m_fname(fname), m_rows(rows), m_cols(cols), m_mode(mode) {
+
+    if (std::filesystem::exists(fname)) {
+        n_frames = std::filesystem::file_size(fname) / (sizeof(sls_detector_header) + rows * cols * bitdepth / 8);
+    } else {
+        n_frames = 0;
+    }
+
+    if (mode == "r") {
+        fp = fopen(m_fname.string().c_str(), "rb");
+    } else {
+        // if file exists, open in read/write mode (without truncating the file)
+        // if file does not exist, open in write mode
+        if (std::filesystem::exists(fname)) {
+            fp = fopen(m_fname.string().c_str(), "r+b");
+        } else {
+            fp = fopen(m_fname.string().c_str(), "wb");
+        }
+    }
+    if (fp == nullptr) {
+        throw std::runtime_error(LOCATION + "Could not open file for writing");
+    }
+}
+
+size_t SubFile::get_part(std::byte *buffer, size_t frame_index) {
+    if (frame_index >= n_frames) {
+        throw std::runtime_error("Frame number out of range");
+    }
+    fseek(fp, sizeof(sls_detector_header) + (sizeof(sls_detector_header) + bytes_per_part()) * frame_index, // NOLINT
+          SEEK_SET);
+    return fread(buffer, this->bytes_per_part(), 1, this->fp);
+}
+size_t SubFile::write_part(std::byte *buffer, sls_detector_header header, size_t frame_index) {
+    if (frame_index > n_frames) {
+        throw std::runtime_error("Frame number out of range");
+    }
+    fseek(fp, static_cast<int64_t>((sizeof(sls_detector_header) + bytes_per_part()) * frame_index), SEEK_SET);
+    auto wc = fwrite(reinterpret_cast<char *>(&header), sizeof(header), 1, fp);
+    wc += fwrite(buffer, bytes_per_part(), 1, fp);
+
+    return wc;
+}
+
+size_t SubFile::frame_number(size_t frame_index) {
+    sls_detector_header h{};
+    fseek(fp, (sizeof(sls_detector_header) + bytes_per_part()) * frame_index, SEEK_SET); // NOLINT
+    size_t const rc = fread(reinterpret_cast<char *>(&h), sizeof(h), 1, fp);
+    if (rc != 1)
+        throw std::runtime_error(LOCATION + "Could not read header from file");
+
+    return h.frameNumber;
+}
+
+SubFile::~SubFile() {
+    if (fp) {
+        fclose(fp);
+    }
+}
+
+} // namespace aare

src/defs.cpp

@@ -0,0 +1,65 @@
+#include "aare/defs.hpp"
+#include <stdexcept>
+#include <string>
+
+namespace aare {
+
+/**
+ * @brief Convert a DetectorType to a string
+ * @param type DetectorType
+ * @return string representation of the DetectorType
+ */
+template <> std::string toString(DetectorType arg) {
+    switch (arg) {
+    case DetectorType::Jungfrau:
+        return "Jungfrau";
+    case DetectorType::Eiger:
+        return "Eiger";
+    case DetectorType::Mythen3:
+        return "Mythen3";
+    case DetectorType::Moench:
+        return "Moench";
+    case DetectorType::ChipTestBoard:
+        return "ChipTestBoard";
+    default:
+        return "Unknown";
+    }
+}
+
+/**
+ * @brief Convert a string to a DetectorType
+ * @param name string representation of the DetectorType
+ * @return DetectorType
+ * @throw runtime_error if the string does not match any DetectorType
+ */
+template <> DetectorType StringTo(const std::string &arg) {
+    if (arg == "Jungfrau")
+        return DetectorType::Jungfrau;
+    if (arg == "Eiger")
+        return DetectorType::Eiger;
+    if (arg == "Mythen3")
+        return DetectorType::Mythen3;
+    if (arg == "Moench")
+        return DetectorType::Moench;
+    if (arg == "ChipTestBoard")
+        return DetectorType::ChipTestBoard;
+    throw std::runtime_error("Could not decode dector from: \"" + arg + "\"");
+}
+
+/**
+ * @brief Convert a string to a TimingMode
+ * @param mode string representation of the TimingMode
+ * @return TimingMode
+ * @throw runtime_error if the string does not match any TimingMode
+ */
+template <> TimingMode StringTo(const std::string &arg) {
+    if (arg == "auto")
+        return TimingMode::Auto;
+    if (arg == "trigger")
+        return TimingMode::Trigger;
+    throw std::runtime_error("Could not decode timing mode from: \"" + arg + "\"");
+}
+
+// template <> TimingMode StringTo<TimingMode>(std::string mode);
+
+} // namespace aare

src/defs.test.cpp

@@ -0,0 +1,42 @@
+#include "aare/defs.hpp"
+// #include "aare/utils/floats.hpp"
+
+#include <catch2/catch_test_macros.hpp>
+#include <string>
+TEST_CASE("Enum to string conversion") {
+    // By the way I don't think the enum string conversions should be in the defs.hpp file
+    // but let's use this to show a test
+    REQUIRE(toString(aare::DetectorType::Jungfrau) == "Jungfrau");
+}
+
+TEST_CASE("Cluster creation") {
+    aare::Cluster c(13, 15);
+    REQUIRE(c.cluster_sizeX == 13);
+    REQUIRE(c.cluster_sizeY == 15);
+    REQUIRE(c.dt == aare::Dtype(typeid(int32_t)));
+    REQUIRE(c.data() != nullptr);
+
+    aare::Cluster c2(c);
+    REQUIRE(c2.cluster_sizeX == 13);
+    REQUIRE(c2.cluster_sizeY == 15);
+    REQUIRE(c2.dt == aare::Dtype(typeid(int32_t)));
+    REQUIRE(c2.data() != nullptr);
+}
+
+// TEST_CASE("cluster set and get data") {
+
+//     aare::Cluster c2(33, 44, aare::Dtype(typeid(double)));
+//     REQUIRE(c2.cluster_sizeX == 33);
+//     REQUIRE(c2.cluster_sizeY == 44);
+//     REQUIRE(c2.dt == aare::Dtype::DOUBLE);
+//     double v = 3.14;
+//     c2.set<double>(0, v);
+//     double v2 = c2.get<double>(0);
+//     REQUIRE(aare::compare_floats<double>(v, v2));
+
+//     c2.set<double>(33 * 44 - 1, 123.11);
+//     double v3 = c2.get<double>(33 * 44 - 1);
+//     REQUIRE(aare::compare_floats<double>(123.11, v3));
+
+//     REQUIRE_THROWS_AS(c2.set(0, 1), std::invalid_argument); // set int to double
+// }

tests/CMakeLists.txt

@@ -0,0 +1,44 @@
+
+if (AARE_FETCH_CATCH)
+    FetchContent_Declare(
+        Catch2
+        GIT_SHALLOW    TRUE
+        GIT_REPOSITORY https://github.com/catchorg/Catch2.git
+        GIT_TAG        v3.5.3
+    )
+    FetchContent_MakeAvailable(Catch2)
+else()
+    find_package(Catch2 3 REQUIRED)
+endif()
+
+
+list(APPEND CMAKE_MODULE_PATH ${Catch2_SOURCE_DIR}/extras)
+
+add_executable(tests test.cpp)
+target_link_libraries(tests PRIVATE Catch2::Catch2WithMain)
+
+set_target_properties(tests PROPERTIES
+        RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}
+        OUTPUT_NAME run_tests
+)
+
+include(CTest)
+include(Catch)
+catch_discover_tests(tests)
+
+set(TestSources
+        ${CMAKE_CURRENT_SOURCE_DIR}/test.cpp
+    )
+target_sources(tests PRIVATE ${TestSources} )
+
+#Work around to remove, this is not the way to do it =) 
+# target_include_directories(tests PRIVATE ${CMAKE_SOURCE_DIR}/include/common)
+target_link_libraries(tests PRIVATE aare_core aare_compiler_flags)
+
+
+
+#configure a header to pass test file paths
+get_filename_component(TEST_FILE_PATH ${PROJECT_SOURCE_DIR}/data ABSOLUTE)
+configure_file(test_config.hpp.in test_config.hpp)
+target_include_directories(tests PRIVATE ${CMAKE_CURRENT_BINARY_DIR})
+

tests/test.cpp

@@ -0,0 +1,21 @@
+#include "test_config.hpp"
+#include <catch2/catch_test_macros.hpp>
+#include <climits>
+#include <filesystem>
+#include <fstream>
+
+TEST_CASE("Test suite can find data assets") {
+    auto fpath = test_data_path() / "numpy" / "test_numpy_file.npy";
+    REQUIRE(std::filesystem::exists(fpath));
+}
+
+TEST_CASE("Test suite can open data assets") {
+    auto fpath = test_data_path() / "numpy" / "test_numpy_file.npy";
+    auto f = std::ifstream(fpath, std::ios::binary);
+    REQUIRE(f.is_open());
+}
+
+TEST_CASE("Test float32 and char8") {
+    REQUIRE(sizeof(float) == 4);
+    REQUIRE(CHAR_BIT == 8);
+}

tests/test_config.hpp.in

@@ -0,0 +1,12 @@
+#pragma once
+#include <filesystem>
+#include <cstdlib>
+
+
+inline auto test_data_path(){
+    if(const char* env_p = std::getenv("AARE_TEST_DATA")){
+        return std::filesystem::path(env_p);
+    }else{
+        throw std::runtime_error("AARE_TEST_DATA_PATH not set");
+    }
+}