Multi threaded fitting and returning chi2 (#132)

Co-authored-by: Patrick <patrick.sieberer@psi.ch>
Co-authored-by: JulianHeymes <julian.heymes@psi.ch>
Co-authored-by: Dhanya Thattil <dhanya.thattil@psi.ch>

.clang-format

@@ -0,0 +1,7 @@
+BasedOnStyle: LLVM
+IndentWidth: 4
+
+UseTab: Never
+ColumnLimit: 80
+AlignConsecutiveAssignments: false 
+AlignConsecutiveMacros: true

.github/workflows/build_and_deploy_conda.yml

@@ -0,0 +1,42 @@
+name: Build pkgs and deploy if on main
+
+on:
+  push:
+    branches:
+     - main
+
+jobs:
+  build:
+    strategy:
+      fail-fast: false
+      matrix:
+        platform: [ubuntu-latest, ] # macos-12, windows-2019]
+        python-version: ["3.12",]
+
+    runs-on: ${{ matrix.platform }}
+
+    # The setup-miniconda action needs this to activate miniconda
+    defaults:
+      run:
+        shell: "bash -l {0}"
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Get conda
+        uses: conda-incubator/setup-miniconda@v3.0.4
+        with:
+          python-version: ${{ matrix.python-version }}
+          channels: conda-forge
+
+      - name: Prepare
+        run: conda install conda-build=24.9 conda-verify pytest anaconda-client
+
+      - name: Enable upload
+        run: conda config --set anaconda_upload yes
+
+      - name: Build
+        env:
+          CONDA_TOKEN: ${{ secrets.CONDA_TOKEN }}
+        run: conda build conda-recipe --user slsdetectorgroup --token ${CONDA_TOKEN}
+

.github/workflows/build_conda.yml

@@ -0,0 +1,40 @@
+name: Build pkgs and deploy if on main
+
+on:
+  push:
+    branches:
+     - developer
+
+jobs:
+  build:
+    strategy:
+      fail-fast: false
+      matrix:
+        platform: [ubuntu-latest, ] # macos-12, windows-2019]
+        python-version: ["3.12",]
+
+    runs-on: ${{ matrix.platform }}
+
+    # The setup-miniconda action needs this to activate miniconda
+    defaults:
+      run:
+        shell: "bash -l {0}"
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Get conda
+        uses: conda-incubator/setup-miniconda@v3.0.4
+        with:
+          python-version: ${{ matrix.python-version }}
+          channels: conda-forge
+
+      - name: Prepare
+        run: conda install conda-build=24.9 conda-verify pytest anaconda-client
+
+      - name: Disable upload
+        run: conda config --set anaconda_upload no
+
+      - name: Build
+        run: conda build conda-recipe 
+

.github/workflows/build_docs.yml

@@ -0,0 +1,68 @@
+name: Build the package using cmake then documentation
+
+on:
+  workflow_dispatch:
+  push:
+    
+
+
+permissions:
+  contents: read
+  pages: write
+  id-token: write
+
+jobs:
+  build:
+    strategy:
+      fail-fast: false
+      matrix:
+        platform: [ubuntu-latest, ] # macos-12, windows-2019]
+        python-version: ["3.12",]
+
+    runs-on: ${{ matrix.platform }}
+
+    # The setup-miniconda action needs this to activate miniconda
+    defaults:
+      run:
+        shell: "bash -l {0}"
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Get conda
+        uses: conda-incubator/setup-miniconda@v3.0.4
+        with:
+          python-version: ${{ matrix.python-version }}
+          channels: conda-forge
+
+      - name: Prepare
+        run: conda install doxygen sphinx=7.1.2 breathe pybind11 sphinx_rtd_theme furo nlohmann_json zeromq fmt numpy
+
+      - name: Build library 
+        run: |
+            mkdir build 
+            cd build 
+            cmake .. -DAARE_SYSTEM_LIBRARIES=ON -DAARE_DOCS=ON 
+            make -j 2 
+            make docs
+
+      - name: Upload static files as artifact
+        id: deployment
+        uses: actions/upload-pages-artifact@v3 
+        with:
+         path: build/docs/html/
+  deploy:
+    environment:
+      name: github-pages
+      url: ${{ steps.deployment.outputs.page_url }}
+    runs-on: ubuntu-latest
+    needs: build
+    if: github.ref == 'refs/heads/main'
+    steps:
+      - name: Deploy to GitHub Pages
+        id: deployment
+        uses: actions/deploy-pages@v4
+
+
+
+

.gitignore

@@ -0,0 +1,24 @@
+install/
+.cproject
+.project
+bin/
+.settings
+*.aux
+*.log
+*.out
+*.toc
+*.o
+*.so
+.*
+build/
+RELEASE.txt
+Testing/
+
+ctbDict.cpp
+ctbDict.h
+
+
+
+*.pyc
+*/__pycache__/*
+

CMakeLists.txt

@@ -0,0 +1,476 @@
+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_std_fs_lib()
+message(STATUS "Extra linking to fs lib:${STD_FS_LIB}")
+
+set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH})
+
+
+# General options
+option(AARE_PYTHON_BINDINGS "Build python bindings" ON)
+option(AARE_TESTS "Build tests" OFF)
+option(AARE_BENCHMARKS "Build benchmarks" OFF)
+option(AARE_EXAMPLES "Build examples" OFF)
+option(AARE_IN_GITHUB_ACTIONS "Running in Github Actions" OFF)
+option(AARE_DOCS "Build documentation" OFF)
+option(AARE_VERBOSE "Verbose output" OFF)
+option(AARE_CUSTOM_ASSERT "Use custom assert" OFF)
+option(AARE_INSTALL_PYTHONEXT "Install the python extension in the install tree under CMAKE_INSTALL_PREFIX/aare/" OFF)
+option(AARE_ASAN "Enable AddressSanitizer" OFF)
+
+# Configure which of the dependencies to use FetchContent for
+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(AARE_FETCH_LMFIT "Use FetchContent to download lmfit" ON)
+
+
+#Convenience option to use system libraries only (no FetchContent)
+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()
+
+if(AARE_VERBOSE)
+    add_compile_definitions(AARE_VERBOSE)
+endif()
+
+if(AARE_CUSTOM_ASSERT)
+    add_compile_definitions(AARE_CUSTOM_ASSERT)
+endif()
+
+if(AARE_BENCHMARKS)
+    add_subdirectory(benchmarks)
+endif()
+
+
+set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
+
+if(AARE_FETCH_LMFIT)
+    set(lmfit_patch git apply ${CMAKE_CURRENT_SOURCE_DIR}/patches/lmfit.patch)
+    FetchContent_Declare(
+        lmfit
+        GIT_REPOSITORY https://jugit.fz-juelich.de/mlz/lmfit.git
+        GIT_TAG        main
+        PATCH_COMMAND ${lmfit_patch}
+        UPDATE_DISCONNECTED 1
+        EXCLUDE_FROM_ALL 1
+    )
+    #Disable what we don't need from lmfit
+    set(BUILD_TESTING OFF CACHE BOOL "")
+    set(LMFIT_CPPTEST OFF CACHE BOOL "")
+    set(LIB_MAN OFF CACHE BOOL "")
+    set(LMFIT_CPPTEST OFF CACHE BOOL "")
+    set(BUILD_SHARED_LIBS OFF CACHE BOOL "")
+
+    
+    FetchContent_MakeAvailable(lmfit)
+    set_property(TARGET lmfit PROPERTY POSITION_INDEPENDENT_CODE ON)
+else()
+    find_package(lmfit REQUIRED)
+endif()
+
+
+if(AARE_FETCH_ZMQ)
+    # Fetchcontent_Declare is deprecated need to find a way to update this
+    # for now setting the policy to old is enough
+    if (${CMAKE_VERSION} VERSION_GREATER_EQUAL "3.30")
+        cmake_policy(SET CMP0169 OLD)
+    endif()
+    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
+    )
+    set(FMT_INSTALL ON CACHE BOOL "")
+    # set(FMT_CMAKE_DIR "")
+    FetchContent_MakeAvailable(fmt)
+    set_property(TARGET fmt PROPERTY POSITION_INDEPENDENT_CODE ON)
+    install(TARGETS fmt 
+        EXPORT ${project}-targets
+        LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
+        ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
+        RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}
+        INCLUDES DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}   
+    )  
+else()
+      find_package(fmt 6 REQUIRED)
+endif()
+
+
+if (AARE_FETCH_JSON)
+    FetchContent_Declare(
+        json 
+        URL https://github.com/nlohmann/json/releases/download/v3.11.3/json.tar.xz
+    )
+    set(JSON_Install ON CACHE BOOL "")
+    FetchContent_MakeAvailable(json)
+    set(NLOHMANN_JSON_TARGET_NAME nlohmann_json)
+    
+    install(
+        TARGETS nlohmann_json
+        EXPORT "${TARGETS_EXPORT_NAME}"
+    )
+    message(STATUS "target: ${NLOHMANN_JSON_TARGET_NAME}")
+else()
+    find_package(nlohmann_json 3.11.3 REQUIRED)
+endif()
+
+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")
+endif()
+
+# Common flags for GCC and Clang
+target_compile_options(
+    aare_compiler_flags 
+    INTERFACE 
+        -Wall 
+        -Wextra 
+        -pedantic 
+        -Wshadow 
+        -Wold-style-cast
+        -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() #GCC/Clang specific
+
+
+if(AARE_ASAN)
+    message(STATUS "AddressSanitizer enabled")
+    target_compile_options(
+        aare_compiler_flags 
+        INTERFACE 
+            -fsanitize=address,undefined,pointer-compare 
+            -fno-omit-frame-pointer
+    )
+    target_link_libraries(
+        aare_compiler_flags 
+        INTERFACE 
+            -fsanitize=address,undefined,pointer-compare 
+            -fno-omit-frame-pointer
+    )
+endif()
+
+
+
+
+
+if(AARE_TESTS)
+    enable_testing()
+    add_subdirectory(tests)
+endif()
+
+###------------------------------------------------------------------------------MAIN LIBRARY
+###------------------------------------------------------------------------------------------
+
+set(PUBLICHEADERS
+    include/aare/ArrayExpr.hpp
+    include/aare/ClusterFinder.hpp
+    include/aare/ClusterFile.hpp
+    include/aare/CtbRawFile.hpp
+    include/aare/ClusterVector.hpp
+    include/aare/decode.hpp
+    include/aare/defs.hpp
+    include/aare/Dtype.hpp
+    include/aare/File.hpp
+    include/aare/Fit.hpp
+    include/aare/FileInterface.hpp
+    include/aare/Frame.hpp
+    include/aare/geo_helpers.hpp
+    include/aare/NDArray.hpp
+    include/aare/NDView.hpp
+    include/aare/NumpyFile.hpp
+    include/aare/NumpyHelpers.hpp
+    include/aare/Pedestal.hpp
+    include/aare/PixelMap.hpp
+    include/aare/RawFile.hpp
+    include/aare/RawMasterFile.hpp
+    include/aare/RawSubFile.hpp
+    include/aare/VarClusterFinder.hpp
+    include/aare/utils/task.hpp
+
+)
+
+
+set(SourceFiles
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/CtbRawFile.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFile.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/defs.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/decode.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/File.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/Fit.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/geo_helpers.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/PixelMap.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/RawSubFile.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/utils/task.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 
+    nlohmann_json::nlohmann_json 
+    ${STD_FS_LIB} # from helpers.cmake
+    PRIVATE 
+    aare_compiler_flags 
+    "$<BUILD_INTERFACE:lmfit>"
+
+)
+
+set_target_properties(aare_core PROPERTIES
+    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/geo_helpers.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.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/ClusterVector.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
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/utils/task.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 aare-core aare-compiler-flags )
+    include(cmake/package_config.cmake)
+endif()

README.md

@@ -1,2 +1,71 @@
 # aare
 Data analysis library for PSI hybrid detectors
+
+
+## Build and install
+
+Prerequisites
+- cmake >= 3.14
+- C++17 compiler (gcc >= 8)
+- python >= 3.10
+
+### Development install (for Python)
+
+```bash
+git clone git@github.com:slsdetectorgroup/aare.git --branch=v1 #or using http...
+mkdir build
+cd build
+
+#configure using cmake
+cmake ../aare
+
+#build (replace 4 with the number of threads you want to use)
+make -j4 
+```
+
+Now you can use the Python module from your build directory
+
+```python
+import aare
+f = aare.File('Some/File/I/Want_to_open_master_0.json')
+```
+
+To run form other folders either add the path to your conda environment using conda-build or add it to your PYTHONPATH
+
+
+### Install using conda/mamba
+
+```bash
+#enable your env first!
+conda install aare=2024.10.29.dev0 -c slsdetectorgroup
+```
+
+### Install to a custom location and use in your project
+
+Working example in: https://github.com/slsdetectorgroup/aare-examples
+
+```bash
+#build and install aare 
+git clone git@github.com:slsdetectorgroup/aare.git --branch=v1 #or using http...
+mkdir build
+cd build
+
+#configure using cmake
+cmake ../aare -DCMAKE_INSTALL_PREFIX=/where/to/put/aare
+
+#build (replace 4 with the number of threads you want to use)
+make -j4 
+
+#install
+make install
+
+
+#Now configure your project
+ cmake .. -DCMAKE_PREFIX_PATH=SOME_PATH
+```
+
+### Local build of conda pkgs
+
+```bash
+conda build . --variants="{python: [3.11, 3.12, 3.13]}"
+```

benchmarks/CMakeLists.txt

@@ -0,0 +1,11 @@
+find_package(benchmark REQUIRED)
+
+add_executable(ndarray_benchmark ndarray_benchmark.cpp)
+
+target_link_libraries(ndarray_benchmark benchmark::benchmark aare_core aare_compiler_flags)
+# target_link_libraries(tests PRIVATE aare_core aare_compiler_flags)
+
+set_target_properties(ndarray_benchmark PROPERTIES
+        RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}
+        # OUTPUT_NAME run_tests
+)

benchmarks/ndarray_benchmark.cpp

@@ -0,0 +1,136 @@
+#include <benchmark/benchmark.h>
+#include "aare/NDArray.hpp"
+
+
+using aare::NDArray;
+
+constexpr ssize_t size = 1024;
+class TwoArrays : public benchmark::Fixture {
+public:
+    NDArray<int,2> a{{size,size},0};
+    NDArray<int,2> b{{size,size},0};
+  void SetUp(::benchmark::State& state) {
+    for(uint32_t i = 0; i < size; i++){
+      for(uint32_t j = 0; j < size; j++){
+        a(i, j)= i*j+1;
+        b(i, j)= i*j+1;
+      }
+    }
+  }
+
+  // void TearDown(::benchmark::State& state) {
+  // }
+};
+
+
+
+
+BENCHMARK_F(TwoArrays, AddWithOperator)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res = a+b;
+    benchmark::DoNotOptimize(res);
+  }
+}
+BENCHMARK_F(TwoArrays, AddWithIndex)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res(a.shape());
+    for (uint32_t i = 0; i < a.size(); i++) {
+        res(i) = a(i) + b(i);
+    }
+    benchmark::DoNotOptimize(res);
+  }
+}
+
+BENCHMARK_F(TwoArrays, SubtractWithOperator)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res = a-b;
+    benchmark::DoNotOptimize(res);
+  }
+}
+BENCHMARK_F(TwoArrays, SubtractWithIndex)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res(a.shape());
+    for (uint32_t i = 0; i < a.size(); i++) {
+        res(i) = a(i) - b(i);
+    }
+    benchmark::DoNotOptimize(res);
+  }
+}
+
+BENCHMARK_F(TwoArrays, MultiplyWithOperator)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res = a*b;
+    benchmark::DoNotOptimize(res);
+  }
+}
+BENCHMARK_F(TwoArrays, MultiplyWithIndex)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res(a.shape());
+    for (uint32_t i = 0; i < a.size(); i++) {
+        res(i) = a(i) * b(i);
+    }
+    benchmark::DoNotOptimize(res);
+  }
+}
+
+BENCHMARK_F(TwoArrays, DivideWithOperator)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res = a/b;
+    benchmark::DoNotOptimize(res);
+  }
+}
+BENCHMARK_F(TwoArrays, DivideWithIndex)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res(a.shape());
+    for (uint32_t i = 0; i < a.size(); i++) {
+        res(i) = a(i) / b(i);
+    }
+    benchmark::DoNotOptimize(res);
+  }
+}
+
+BENCHMARK_F(TwoArrays, FourAddWithOperator)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res = a+b+a+b;
+    benchmark::DoNotOptimize(res);
+  }
+}
+BENCHMARK_F(TwoArrays, FourAddWithIndex)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res(a.shape());
+    for (uint32_t i = 0; i < a.size(); i++) {
+        res(i) = a(i) + b(i) + a(i) + b(i);
+    }
+    benchmark::DoNotOptimize(res);
+  }
+}
+
+BENCHMARK_F(TwoArrays, MultiplyAddDivideWithOperator)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res = a*a+b/a;
+    benchmark::DoNotOptimize(res);
+  }
+}
+BENCHMARK_F(TwoArrays, MultiplyAddDivideWithIndex)(benchmark::State& st) {
+   for (auto _ : st) {
+    // This code gets timed
+    NDArray<int,2> res(a.shape());
+    for (uint32_t i = 0; i < a.size(); i++) {
+        res(i) = a(i) * a(i) + b(i) / a(i);
+    }
+    benchmark::DoNotOptimize(res);
+  }
+}
+
+BENCHMARK_MAIN();

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,46 @@
+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()
+
+function(set_std_fs_lib)
+# from pybind11
+# Check if we need to add -lstdc++fs or -lc++fs or nothing
+if(DEFINED CMAKE_CXX_STANDARD AND CMAKE_CXX_STANDARD LESS 17)
+  set(STD_FS_NO_LIB_NEEDED TRUE)
+elseif(MSVC)
+  set(STD_FS_NO_LIB_NEEDED TRUE)
+else()
+  file(
+    WRITE ${CMAKE_CURRENT_BINARY_DIR}/main.cpp
+    "#include <filesystem>\nint main(int argc, char ** argv) {\n  std::filesystem::path p(argv[0]);\n  return p.string().length();\n}"
+  )
+  try_compile(
+    STD_FS_NO_LIB_NEEDED ${CMAKE_CURRENT_BINARY_DIR}
+    SOURCES ${CMAKE_CURRENT_BINARY_DIR}/main.cpp
+    COMPILE_DEFINITIONS -std=c++17)
+  try_compile(
+    STD_FS_NEEDS_STDCXXFS ${CMAKE_CURRENT_BINARY_DIR}
+    SOURCES ${CMAKE_CURRENT_BINARY_DIR}/main.cpp
+    COMPILE_DEFINITIONS -std=c++17
+    LINK_LIBRARIES stdc++fs)
+  try_compile(
+    STD_FS_NEEDS_CXXFS ${CMAKE_CURRENT_BINARY_DIR}
+    SOURCES ${CMAKE_CURRENT_BINARY_DIR}/main.cpp
+    COMPILE_DEFINITIONS -std=c++17
+    LINK_LIBRARIES c++fs)
+endif()
+
+if(${STD_FS_NEEDS_STDCXXFS})
+  set(STD_FS_LIB stdc++fs PARENT_SCOPE)
+elseif(${STD_FS_NEEDS_CXXFS})
+  set(STD_FS_LIB c++fs PARENT_SCOPE)
+elseif(${STD_FS_NO_LIB_NEEDED})
+  set(STD_FS_LIB "" PARENT_SCOPE)
+else()
+  message(WARNING "Unknown C++17 compiler - not passing -lstdc++fs")
+  set(STD_FS_LIB "")
+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,28 @@
+# 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@")
+
+# List dependencies
+find_dependency(Threads)
+find_dependency(fmt)
+find_dependency(nlohmann_json)
+
+# 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/conda_build_config.yaml

@@ -0,0 +1,28 @@
+python:
+    - 3.11
+    - 3.11
+    - 3.11
+    - 3.12
+    - 3.12
+    - 3.12
+    - 3.13
+
+
+    
+numpy: 
+    - 1.26
+    - 2.0
+    - 2.1
+    - 1.26
+    - 2.0
+    - 2.1
+    - 2.1
+
+
+zip_keys:
+    - python
+    - numpy
+
+pin_run_as_build:
+  numpy: x.x
+  python: x.x

conda-recipe/meta.yaml

@@ -0,0 +1,53 @@
+package:
+  name: aare
+  version: 2025.2.18 #TODO! how to not duplicate this?
+
+
+
+source:
+  path: ..
+
+build:
+  number: 0
+  script:
+    - unset CMAKE_GENERATOR && {{ PYTHON }} -m pip install . -vv  # [not win]
+    - {{ PYTHON }} -m pip install . -vv  # [win]
+
+requirements:
+  build:
+    - python {{python}}
+    - numpy {{ numpy }}
+    - {{ compiler('cxx') }}
+
+
+  host:
+    - cmake
+    - ninja
+    - python {{python}}
+    - numpy {{ numpy }}
+    - pip
+    - scikit-build-core
+    - pybind11 >=2.13.0
+    - fmt
+    - zeromq
+    - nlohmann_json
+    - catch2
+
+  run:
+    - python {{python}}
+    - numpy {{ numpy }}
+
+
+test:
+  imports:
+    - aare
+  # requires:
+  #   - pytest
+  # source_files:
+  #   - tests
+  # commands:
+  #   - pytest tests
+
+about:
+  summary: An example project built with pybind11 and scikit-build.
+  # license_file: LICENSE

docs/CMakeLists.txt

@@ -0,0 +1,56 @@
+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})
+
+
+file(GLOB SPHINX_SOURCE_FILES CONFIGURE_DEPENDS "src/*.rst")
+
+
+
+foreach(filename ${SPHINX_SOURCE_FILES})
+    get_filename_component(fname ${filename} NAME)
+    message(STATUS "Copying ${filename} to ${SPHINX_BUILD}/src/${fname}")
+    configure_file(${filename} "${SPHINX_BUILD}/src/${fname}")  
+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,63 @@
+# 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('..'))
+
+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.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/ClusterFile.rst

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

docs/src/ClusterFinder.rst

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

docs/src/ClusterFinderMT.rst

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

docs/src/ClusterVector.rst

@@ -0,0 +1,6 @@
+ClusterVector
+=============
+
+.. doxygenclass:: aare::ClusterVector
+   :members:
+   :undoc-members:

docs/src/Consume.rst

@@ -0,0 +1,19 @@
+Use from C++
+========================
+
+There are a few different way to use aare in your C++ project. Which one you choose 
+depends on how you intend to work with the library and how you manage your dependencies.
+
+
+Install and use cmake with find_package(aare)
+-------------------------------------------------
+
+https://github.com/slsdetectorgroup/aare-examples
+
+.. include:: _install.rst
+
+
+Use as a submodule
+------------------- 
+
+Coming soon...

docs/src/Dtype.rst

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

docs/src/File.rst

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

docs/src/Frame.rst

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

docs/src/Installation.rst

@@ -0,0 +1,106 @@
+****************
+Installation
+****************
+
+.. attention ::
+
+    - https://cliutils.gitlab.io/modern-cmake/README.html
+
+conda/mamaba
+~~~~~~~~~~~~~~~~~~
+
+This is the recommended way to install aare. Using a package manager makes it easy to 
+switch between versions and is (one of) the most convenient way to install up to date
+dependencies on older distributions.
+
+.. note ::
+
+    aare is developing rapidly. Check for the latest release by 
+    using: **conda search aare -c slsdetectorgroup**
+
+
+.. code-block:: bash
+
+    # Install a specific version: 
+    conda install aare=2024.11.11.dev0 -c slsdetectorgroup
+
+
+cmake build (development install)
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+If you are working on aare or want to test our a version that doesn't yet have
+a conda package. Build using cmake and then run from the build folder.
+
+.. code-block:: bash
+
+    git clone git@github.com:slsdetectorgroup/aare.git --branch=v1 #or using http...
+    mkdir build
+    cd build
+
+    #configure using cmake
+    cmake ../aare
+
+    #build (replace 4 with the number of threads you want to use)
+    make -j4 
+
+
+    # add the build folder to your PYTHONPATH and then you should be able to
+    # import aare in python
+
+cmake build + install and use in your C++ project
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. warning ::
+
+    When building aare with default settings we also include fmt and nlohmann_json.
+    Installation to a custom location is highly recommended.
+
+
+.. note ::
+
+    It is also possible to install aare with conda and then use in your C++ project.
+
+.. include:: _install.rst
+
+
+cmake options
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+For detailed options see the CMakeLists.txt file in the root directory of the project.
+
+.. code-block:: bash
+
+    # usage (or edit with ccmake .)
+    cmake ../aare -DOPTION1=ON -DOPTION2=OFF
+
+
+**AARE_SYSTEM_LIBRARIES "Use system libraries" OFF**
+
+Use system libraries instead of using FetchContent to pull in dependencies. Default option is off.
+
+
+**AARE_PYTHON_BINDINGS "Build python bindings" ON**
+
+Build the Python bindings. Default option is on. 
+
+.. warning ::
+
+    If you have a newer system Python compared to the one in your virtual environment,
+    you might have to pass -DPython_FIND_VIRTUALENV=ONLY to cmake.
+
+**AARE_TESTS "Build tests" OFF**
+
+Build unit tests. Default option is off.
+
+**AARE_EXAMPLES "Build examples" OFF**
+
+**AARE_DOCS "Build documentation" OFF**
+
+Build documentation. Needs doxygen, sphinx and breathe. Default option is off.
+Requires a separate make docs.
+
+**AARE_VERBOSE "Verbose output" OFF**
+
+**AARE_CUSTOM_ASSERT "Use custom assert" OFF**
+
+Enable custom assert macro to check for errors. Default option is off.

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,8 @@
+Pedestal
+=============
+
+
+.. doxygenclass:: aare::Pedestal
+   :members:
+   :undoc-members:
+   :private-members:

docs/src/RawFile.rst

@@ -0,0 +1,8 @@
+RawFile
+===============
+
+
+.. doxygenclass:: aare::RawFile
+   :members:
+   :undoc-members:
+   :private-members:

docs/src/RawMasterFile.rst

@@ -0,0 +1,14 @@
+RawMasterFile
+===============
+
+
+.. doxygenclass:: aare::RawMasterFile
+   :members:
+   :undoc-members:
+   :private-members:
+
+
+.. doxygenclass:: aare::RawFileNameComponents
+   :members:
+   :undoc-members:
+   :private-members:

docs/src/RawSubFile.rst

@@ -0,0 +1,8 @@
+RawSubFile
+===============
+
+
+.. doxygenclass:: aare::RawSubFile
+   :members:
+   :undoc-members:
+   :private-members:

docs/src/Requirements.rst

@@ -0,0 +1,23 @@
+Requirements
+==============================================
+
+- C++17 compiler (gcc 8/clang 7)
+- CMake 3.14+
+
+**Internally used libraries**
+
+.. note ::
+
+    These can also be picked up from the system/conda environment by specifying:
+    -DAARE_SYSTEM_LIBRARIES=ON during the cmake configuration.
+
+- pybind11
+- fmt
+- nlohmann_json
+- ZeroMQ
+
+**Extra dependencies for building documentation**
+
+- Sphinx
+- Breathe
+- Doxygen

docs/src/VarClusterFinder.rst

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

docs/src/_install.rst

@@ -0,0 +1,23 @@
+.. code-block:: bash
+
+    #build and install aare 
+    git clone git@github.com:slsdetectorgroup/aare.git --branch=developer #or using http...
+    mkdir build
+    cd build
+
+    #configure using cmake
+    cmake ../aare -DCMAKE_INSTALL_PREFIX=/where/to/put/aare
+
+    #build (replace 4 with the number of threads you want to use)
+    make -j4 
+
+    #install
+    make install
+
+    #Go to your project
+    cd /your/project/source
+
+    #Now configure your project
+    mkdir build
+    cd build
+    cmake .. -DCMAKE_PREFIX_PATH=SOME_PATH

docs/src/index.rst

@@ -0,0 +1,62 @@
+AARE
+==============================================
+
+.. note :: 
+    
+    **Examples:**
+
+    - `jupyter notebooks <https://github.com/slsdetectorgroup/aare-notebooks>`_
+    - `cmake+install <https://github.com/slsdetectorgroup/aare-examples>`_
+    - `git submodule <https://github.com/slsdetectorgroup/aare-submodule>`_
+
+
+
+
+.. toctree::
+    :caption: Installation
+    :maxdepth: 3
+    
+    Installation
+    Requirements
+    Consume
+
+
+
+
+.. toctree::
+    :caption: Python API
+    :maxdepth: 1
+    
+    pyFile
+    pyCtbRawFile
+    pyClusterFile
+    pyClusterVector
+    pyRawFile
+    pyRawMasterFile
+    pyVarClusterFinder
+
+    pyFit
+
+
+.. toctree::
+    :caption: C++ API
+    :maxdepth: 1
+    
+    NDArray
+    NDView
+    Frame
+    File
+    Dtype
+    ClusterFinder
+    ClusterFinderMT
+    ClusterFile
+    ClusterVector
+    Pedestal
+    RawFile
+    RawSubFile
+    RawMasterFile
+    VarClusterFinder
+
+
+
+

docs/src/pyClusterFile.rst

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

docs/src/pyClusterVector.rst

@@ -0,0 +1,33 @@
+ClusterVector
+================
+
+The ClusterVector, holds clusters from the ClusterFinder. Since it is templated
+in C++  we use a suffix indicating the data type in python. The suffix is
+``_i`` for integer, ``_f`` for float, and ``_d`` for double.
+
+At the moment the functionality from python is limited and it is not supported
+to push_back clusters to the vector. The intended use case is to pass it to 
+C++ functions that support the ClusterVector or to view it as a numpy array.
+
+**View ClusterVector as numpy array**
+
+.. code:: python
+
+    from aare import ClusterFile
+    with ClusterFile("path/to/file") as f:
+        cluster_vector = f.read_frame()
+
+    # Create a copy of the cluster data in a numpy array
+    clusters = np.array(cluster_vector)
+
+    # Avoid copying the data by passing copy=False
+    clusters = np.array(cluster_vector, copy = False)
+
+
+.. py:currentmodule:: aare
+
+.. autoclass:: ClusterVector_i
+    :members:
+    :undoc-members:
+    :show-inheritance:
+    :inherited-members:

docs/src/pyCtbRawFile.rst

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

docs/src/pyFile.rst

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

docs/src/pyFit.rst

@@ -0,0 +1,19 @@
+
+Fit
+========
+
+.. py:currentmodule:: aare
+
+
+**Functions**
+
+.. autofunction:: gaus
+
+.. autofunction:: pol1
+
+
+**Fitting**
+
+.. autofunction:: fit_gaus
+
+.. autofunction:: fit_pol1

docs/src/pyRawFile.rst

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

docs/src/pyRawMasterFile.rst

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

docs/src/pyVarClusterFinder.rst

@@ -0,0 +1,10 @@
+VarClusterFinder
+===================
+
+.. py:currentmodule:: aare
+
+.. autoclass:: VarClusterFinder
+    :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/ArrayExpr.hpp

@@ -0,0 +1,99 @@
+#pragma once
+#include <cstdint> //int64_t
+#include <cstddef> //size_t
+#include <array> 
+
+#include <cassert>
+namespace aare {
+
+template <typename E, int64_t Ndim> class ArrayExpr {
+  public:
+    static constexpr bool is_leaf = false;
+
+    auto operator[](size_t i) const { return static_cast<E const &>(*this)[i]; }
+    auto operator()(size_t i) const { return static_cast<E const &>(*this)[i]; }
+    auto size() const { return static_cast<E const &>(*this).size(); }
+    std::array<int64_t, Ndim> shape() const { return static_cast<E const &>(*this).shape(); }
+};
+
+template <typename A, typename B, int64_t Ndim>
+class ArrayAdd : public ArrayExpr<ArrayAdd<A, B, Ndim>, Ndim> {
+    const A &arr1_;
+    const B &arr2_;
+
+  public:
+    ArrayAdd(const A &arr1, const B &arr2) : arr1_(arr1), arr2_(arr2) {
+        assert(arr1.size() == arr2.size());
+    }
+    auto operator[](int i) const { return arr1_[i] + arr2_[i]; }
+    size_t size() const { return arr1_.size(); }
+    std::array<int64_t, Ndim> shape() const { return arr1_.shape(); }
+};
+
+template <typename A, typename B, int64_t Ndim>
+class ArraySub : public ArrayExpr<ArraySub<A, B, Ndim>, Ndim> {
+    const A &arr1_;
+    const B &arr2_;
+
+  public:
+    ArraySub(const A &arr1, const B &arr2) : arr1_(arr1), arr2_(arr2) {
+        assert(arr1.size() == arr2.size());
+    }
+    auto operator[](int i) const { return arr1_[i] - arr2_[i]; }
+    size_t size() const { return arr1_.size(); }
+    std::array<int64_t, Ndim> shape() const { return arr1_.shape(); }
+};
+
+template <typename A, typename B, int64_t Ndim>
+class ArrayMul : public ArrayExpr<ArrayMul<A, B, Ndim>,Ndim> {
+    const A &arr1_;
+    const B &arr2_;
+
+  public:
+    ArrayMul(const A &arr1, const B &arr2) : arr1_(arr1), arr2_(arr2) {
+        assert(arr1.size() == arr2.size());
+    }
+    auto operator[](int i) const { return arr1_[i] * arr2_[i]; }
+    size_t size() const { return arr1_.size(); }
+    std::array<int64_t, Ndim> shape() const { return arr1_.shape(); }
+};
+
+template <typename A, typename B, int64_t Ndim>
+class ArrayDiv : public ArrayExpr<ArrayDiv<A, B, Ndim>, Ndim> {
+    const A &arr1_;
+    const B &arr2_;
+
+  public:
+    ArrayDiv(const A &arr1, const B &arr2) : arr1_(arr1), arr2_(arr2) {
+        assert(arr1.size() == arr2.size());
+    }
+    auto operator[](int i) const { return arr1_[i] / arr2_[i]; }
+    size_t size() const { return arr1_.size(); }
+    std::array<int64_t, Ndim> shape() const { return arr1_.shape(); }
+};
+
+
+
+template <typename A, typename B, int64_t Ndim>
+auto operator+(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
+    return ArrayAdd<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
+}
+
+template <typename A, typename B, int64_t Ndim>
+auto operator-(const ArrayExpr<A,Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
+    return ArraySub<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
+}
+
+template <typename A, typename B, int64_t Ndim>
+auto operator*(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
+    return ArrayMul<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
+}
+
+template <typename A, typename B, int64_t Ndim>
+auto operator/(const ArrayExpr<A, Ndim> &arr1, const ArrayExpr<B, Ndim> &arr2) {
+    return ArrayDiv<ArrayExpr<A, Ndim>, ArrayExpr<B, Ndim>, Ndim>(arr1, arr2);
+}
+
+
+
+} // namespace aare

include/aare/CircularFifo.hpp

@@ -0,0 +1,97 @@
+#pragma once
+
+#include <chrono>
+#include <fmt/color.h>
+#include <fmt/format.h>
+#include <memory>
+#include <thread>
+
+#include "aare/ProducerConsumerQueue.hpp"
+
+namespace aare {
+
+template <class ItemType> class CircularFifo {
+    uint32_t fifo_size;
+    aare::ProducerConsumerQueue<ItemType> free_slots;
+    aare::ProducerConsumerQueue<ItemType> filled_slots;
+
+  public:
+    CircularFifo() : CircularFifo(100){};
+    CircularFifo(uint32_t size) : fifo_size(size), free_slots(size + 1), filled_slots(size + 1) {
+
+        // TODO! how do we deal with alignment for writing? alignas???
+        // Do we give the user a chance to provide memory locations?
+        // Templated allocator?
+        for (size_t i = 0; i < fifo_size; ++i) {
+            free_slots.write(ItemType{});
+        }
+    }
+
+    bool next() {
+        // TODO! avoid default constructing ItemType
+        ItemType it;
+        if (!filled_slots.read(it))
+            return false;
+        if (!free_slots.write(std::move(it)))
+            return false;
+        return true;
+    }
+
+    ~CircularFifo() {}
+
+    using value_type = ItemType;
+
+    auto numFilledSlots() const noexcept { return filled_slots.sizeGuess(); }
+    auto numFreeSlots() const noexcept { return free_slots.sizeGuess(); }
+    auto isFull() const noexcept { return filled_slots.isFull(); }
+
+    ItemType pop_free() {
+        ItemType v;
+        while (!free_slots.read(v))
+            ;
+        return std::move(v);
+        // return v;
+    }
+
+    bool try_pop_free(ItemType &v) { return free_slots.read(v); }
+
+    ItemType pop_value(std::chrono::nanoseconds wait, std::atomic<bool> &stopped) {
+        ItemType v;
+        while (!filled_slots.read(v) && !stopped) {
+            std::this_thread::sleep_for(wait);
+        }
+        return std::move(v);
+    }
+
+    ItemType pop_value() {
+        ItemType v;
+        while (!filled_slots.read(v))
+            ;
+        return std::move(v);
+    }
+
+    ItemType *frontPtr() { return filled_slots.frontPtr(); }
+
+    // TODO! Add function to move item from filled to free to be used
+    // with the frontPtr function
+
+    template <class... Args> void push_value(Args &&...recordArgs) {
+        while (!filled_slots.write(std::forward<Args>(recordArgs)...))
+            ;
+    }
+
+    template <class... Args> bool try_push_value(Args &&...recordArgs) {
+        return filled_slots.write(std::forward<Args>(recordArgs)...);
+    }
+
+    template <class... Args> void push_free(Args &&...recordArgs) {
+        while (!free_slots.write(std::forward<Args>(recordArgs)...))
+            ;
+    }
+
+    template <class... Args> bool try_push_free(Args &&...recordArgs) {
+        return free_slots.write(std::forward<Args>(recordArgs)...);
+    }
+};
+
+} // namespace aare

include/aare/ClusterCollector.hpp

@@ -0,0 +1,52 @@
+#pragma once
+#include <atomic>
+#include <thread>
+
+#include "aare/ProducerConsumerQueue.hpp"
+#include "aare/ClusterVector.hpp"
+#include "aare/ClusterFinderMT.hpp"
+
+namespace aare {
+
+class ClusterCollector{
+        ProducerConsumerQueue<ClusterVector<int>>* m_source;
+        std::atomic<bool> m_stop_requested{false};
+        std::atomic<bool> m_stopped{true};
+        std::chrono::milliseconds m_default_wait{1};
+        std::thread m_thread;
+        std::vector<ClusterVector<int>> m_clusters;
+
+    void process(){
+        m_stopped = false;
+        fmt::print("ClusterCollector started\n");
+        while (!m_stop_requested || !m_source->isEmpty()) { 
+            if (ClusterVector<int> *clusters = m_source->frontPtr();
+                clusters != nullptr) {
+                m_clusters.push_back(std::move(*clusters));
+                m_source->popFront();
+            }else{
+                std::this_thread::sleep_for(m_default_wait);
+            }
+        }
+        fmt::print("ClusterCollector stopped\n");
+        m_stopped = true;
+    }
+
+    public:
+        ClusterCollector(ClusterFinderMT<uint16_t, double, int32_t>* source){
+            m_source = source->sink();
+            m_thread = std::thread(&ClusterCollector::process, this);
+        }
+        void stop(){
+            m_stop_requested = true;
+            m_thread.join();
+        }
+        std::vector<ClusterVector<int>> steal_clusters(){
+            if(!m_stopped){
+                throw std::runtime_error("ClusterCollector is still running");
+            }
+            return std::move(m_clusters);
+        }
+};
+
+} // namespace aare

include/aare/ClusterFile.hpp

@@ -0,0 +1,135 @@
+#pragma once
+
+#include "aare/ClusterVector.hpp"
+#include "aare/NDArray.hpp"
+#include "aare/defs.hpp"
+#include <filesystem>
+#include <fstream>
+
+namespace aare {
+
+struct Cluster3x3 {
+    int16_t x;
+    int16_t y;
+    int32_t data[9];
+};
+
+typedef enum {
+    cBottomLeft = 0,
+    cBottomRight = 1,
+    cTopLeft = 2,
+    cTopRight = 3
+} corner;
+
+typedef enum {
+    pBottomLeft = 0,
+    pBottom = 1,
+    pBottomRight = 2,
+    pLeft = 3,
+    pCenter = 4,
+    pRight = 5,
+    pTopLeft = 6,
+    pTop = 7,
+    pTopRight = 8
+} pixel;
+
+struct Eta2 {
+    double x;
+    double y;
+    corner c;
+};
+
+struct ClusterAnalysis {
+    uint32_t c;
+    int32_t tot;
+    double etax;
+    double etay;
+};
+
+/*
+Binary cluster file. Expects data to be layed out as:
+int32_t frame_number
+uint32_t number_of_clusters
+int16_t x, int16_t y, int32_t data[9] x number_of_clusters
+int32_t frame_number
+uint32_t number_of_clusters
+....
+*/
+
+/**
+ * @brief Class to read and write cluster files
+ * Expects data to be laid out as:
+ *
+ *
+ *       int32_t frame_number
+ *       uint32_t number_of_clusters
+ *       int16_t x, int16_t y, int32_t data[9] x number_of_clusters
+ *       int32_t frame_number
+ *       uint32_t number_of_clusters
+ *       etc.
+ */
+class ClusterFile {
+    FILE *fp{};
+    uint32_t m_num_left{};
+    size_t m_chunk_size{};
+    const std::string m_mode;
+
+  public:
+    /**
+     * @brief Construct a new Cluster File object
+     * @param fname path to the file
+     * @param chunk_size number of clusters to read at a time when iterating
+     * over the file
+     * @param mode mode to open the file in. "r" for reading, "w" for writing,
+     * "a" for appending
+     * @throws std::runtime_error if the file could not be opened
+     */
+    ClusterFile(const std::filesystem::path &fname, size_t chunk_size = 1000,
+                const std::string &mode = "r");
+    
+    
+    ~ClusterFile();
+
+    /**
+     * @brief Read n_clusters clusters from the file discarding frame numbers.
+     * If EOF is reached the returned vector will have less than n_clusters
+     * clusters
+     */
+    ClusterVector<int32_t> read_clusters(size_t n_clusters);
+
+    /**
+     * @brief Read a single frame from the file and return the clusters. The
+     * cluster vector will have the frame number set.
+     * @throws std::runtime_error if the file is not opened for reading or the file pointer not
+     * at the beginning of a frame
+     */
+    ClusterVector<int32_t> read_frame();
+
+
+    void write_frame(const ClusterVector<int32_t> &clusters);
+    
+    // Need to be migrated to support NDArray and return a ClusterVector
+    // std::vector<Cluster3x3>
+    // read_cluster_with_cut(size_t n_clusters, double *noise_map, int nx, int ny);
+
+    /**
+     * @brief Return the chunk size
+     */
+    size_t chunk_size() const { return m_chunk_size; }
+    
+    
+    /**
+     * @brief Close the file. If not closed the file will be closed in the destructor
+     */
+    void close();
+};
+
+int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad,
+                 double *eta2x, double *eta2y, double *eta3x, double *eta3y);
+int analyze_cluster(Cluster3x3 &cl, int32_t *t2, int32_t *t3, char *quad,
+                    double *eta2x, double *eta2y, double *eta3x, double *eta3y);
+
+NDArray<double, 2> calculate_eta2(ClusterVector<int> &clusters);
+Eta2 calculate_eta2(Cluster3x3 &cl);
+
+} // namespace aare

include/aare/ClusterFileSink.hpp

@@ -0,0 +1,56 @@
+#pragma once
+#include <atomic>
+#include <filesystem>
+#include <thread>
+
+#include "aare/ProducerConsumerQueue.hpp"
+#include "aare/ClusterVector.hpp"
+#include "aare/ClusterFinderMT.hpp"
+
+namespace aare{
+
+class ClusterFileSink{
+    ProducerConsumerQueue<ClusterVector<int>>* m_source;
+    std::atomic<bool> m_stop_requested{false};
+    std::atomic<bool> m_stopped{true};
+    std::chrono::milliseconds m_default_wait{1};
+    std::thread m_thread;
+    std::ofstream m_file;
+
+
+    void process(){
+        m_stopped = false;
+        fmt::print("ClusterFileSink started\n");
+        while (!m_stop_requested || !m_source->isEmpty()) { 
+            if (ClusterVector<int> *clusters = m_source->frontPtr();
+                clusters != nullptr) {
+                // Write clusters to file
+                int32_t frame_number = clusters->frame_number(); //TODO! Should we store frame number already as int?
+                uint32_t num_clusters = clusters->size();
+                m_file.write(reinterpret_cast<const char*>(&frame_number), sizeof(frame_number));
+                m_file.write(reinterpret_cast<const char*>(&num_clusters), sizeof(num_clusters));
+                m_file.write(reinterpret_cast<const char*>(clusters->data()), clusters->size() * clusters->item_size());
+                m_source->popFront();
+            }else{
+                std::this_thread::sleep_for(m_default_wait);
+            }
+        }
+        fmt::print("ClusterFileSink stopped\n");
+        m_stopped = true;
+    }
+
+    public:
+        ClusterFileSink(ClusterFinderMT<uint16_t, double, int32_t>* source, const std::filesystem::path& fname){
+            m_source = source->sink();
+            m_thread = std::thread(&ClusterFileSink::process, this);
+            m_file.open(fname, std::ios::binary);
+        }
+        void stop(){
+            m_stop_requested = true;
+            m_thread.join();
+            m_file.close();
+        }
+};
+
+
+} // namespace aare

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,149 @@
+#pragma once
+#include "aare/ClusterFile.hpp"
+#include "aare/ClusterVector.hpp"
+#include "aare/Dtype.hpp"
+#include "aare/NDArray.hpp"
+#include "aare/NDView.hpp"
+#include "aare/Pedestal.hpp"
+#include "aare/defs.hpp"
+#include <cstddef>
+
+namespace aare {
+
+template <typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double,
+          typename CT = int32_t>
+class ClusterFinder {
+    Shape<2> m_image_size;
+    const int m_cluster_sizeX;
+    const int m_cluster_sizeY;
+    const PEDESTAL_TYPE m_nSigma;
+    const PEDESTAL_TYPE c2;
+    const PEDESTAL_TYPE c3;
+    Pedestal<PEDESTAL_TYPE> m_pedestal;
+    ClusterVector<CT> m_clusters;
+
+  public:
+    /**
+     * @brief Construct a new ClusterFinder object
+     * @param image_size size of the image
+     * @param cluster_size size of the cluster (x, y)
+     * @param nSigma number of sigma above the pedestal to consider a photon
+     * @param capacity initial capacity of the cluster vector
+     *
+     */
+    ClusterFinder(Shape<2> image_size, Shape<2> cluster_size,
+                  PEDESTAL_TYPE nSigma = 5.0, size_t capacity = 1000000)
+        : m_image_size(image_size), m_cluster_sizeX(cluster_size[0]),
+          m_cluster_sizeY(cluster_size[1]),
+          m_nSigma(nSigma),
+          c2(sqrt((m_cluster_sizeY + 1) / 2 * (m_cluster_sizeX + 1) / 2)),
+          c3(sqrt(m_cluster_sizeX * m_cluster_sizeY)),
+          m_pedestal(image_size[0], image_size[1]),
+          m_clusters(m_cluster_sizeX, m_cluster_sizeY, capacity) {};
+
+    void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {
+        m_pedestal.push(frame);
+    }
+
+    NDArray<PEDESTAL_TYPE, 2> pedestal() { return m_pedestal.mean(); }
+    NDArray<PEDESTAL_TYPE, 2> noise() { return m_pedestal.std(); }
+    void clear_pedestal() { m_pedestal.clear(); }
+
+    /**
+     * @brief Move the clusters from the ClusterVector in the ClusterFinder to a
+     * new ClusterVector and return it.
+     * @param realloc_same_capacity if true the new ClusterVector will have the
+     * same capacity as the old one
+     *
+     */
+    ClusterVector<CT> steal_clusters(bool realloc_same_capacity = false) {
+        ClusterVector<CT> tmp = std::move(m_clusters);
+        if (realloc_same_capacity)
+            m_clusters = ClusterVector<CT>(m_cluster_sizeX, m_cluster_sizeY,
+                                           tmp.capacity());
+        else
+            m_clusters = ClusterVector<CT>(m_cluster_sizeX, m_cluster_sizeY);
+        return tmp;
+    }
+    void find_clusters(NDView<FRAME_TYPE, 2> frame, uint64_t frame_number = 0) {
+        // // TODO! deal with even size clusters
+        // // currently 3,3 -> +/- 1
+        // //  4,4 -> +/- 2
+        int dy = m_cluster_sizeY / 2;
+        int dx = m_cluster_sizeX / 2;
+        m_clusters.set_frame_number(frame_number);
+        std::vector<CT> cluster_data(m_cluster_sizeX * m_cluster_sizeY);
+        for (int iy = 0; iy < frame.shape(0); iy++) {
+            for (int ix = 0; ix < frame.shape(1); ix++) {
+
+                PEDESTAL_TYPE max = std::numeric_limits<FRAME_TYPE>::min();
+                PEDESTAL_TYPE total = 0;
+
+                // What can we short circuit here?
+                PEDESTAL_TYPE rms = m_pedestal.std(iy, ix);
+                PEDESTAL_TYPE value = (frame(iy, ix) - m_pedestal.mean(iy, ix));
+
+                if (value < -m_nSigma * rms)
+                    continue; // NEGATIVE_PEDESTAL go to next pixel
+                              // TODO! No pedestal update???
+
+                for (int ir = -dy; ir < dy + 1; ir++) {
+                    for (int ic = -dx; ic < dx + 1; ic++) {
+                        if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
+                            iy + ir >= 0 && iy + ir < frame.shape(0)) {
+                            PEDESTAL_TYPE val =
+                                frame(iy + ir, ix + ic) -
+                                m_pedestal.mean(iy + ir, ix + ic);
+
+                            total += val;
+                            max = std::max(max, val);
+                        }
+                    }
+                }
+
+                if ((max > m_nSigma * rms)) {
+                    if (value < max)
+                        continue; // Not max go to the next pixel
+                                  // but also no pedestal update
+                } else if (total > c3 * m_nSigma * rms) {
+                    // pass
+                } else {
+                    // m_pedestal.push(iy, ix, frame(iy, ix));   // Safe option
+                    m_pedestal.push_fast(iy, ix, frame(iy, ix)); // Assume we have reached n_samples in the pedestal, slight performance improvement
+                    continue; // It was a pedestal value nothing to store
+                }
+
+                // Store cluster
+                if (value == max) {
+                    // Zero out the cluster data
+                    std::fill(cluster_data.begin(), cluster_data.end(), 0);
+
+                    // Fill the cluster data since we have a photon to store
+                    // It's worth redoing the look since most of the time we
+                    // don't have a photon
+                    int i = 0;
+                    for (int ir = -dy; ir < dy + 1; ir++) {
+                        for (int ic = -dx; ic < dx + 1; ic++) {
+                            if (ix + ic >= 0 && ix + ic < frame.shape(1) &&
+                                iy + ir >= 0 && iy + ir < frame.shape(0)) {
+                                CT tmp =
+                                    static_cast<CT>(frame(iy + ir, ix + ic)) -
+                                    m_pedestal.mean(iy + ir, ix + ic);
+                                cluster_data[i] =
+                                    tmp; // Watch for out of bounds access
+                                i++;
+                            }
+                        }
+                    }
+
+                    // Add the cluster to the output ClusterVector
+                    m_clusters.push_back(
+                        ix, iy,
+                        reinterpret_cast<std::byte *>(cluster_data.data()));
+                }
+            }
+        }
+    }
+};
+
+} // namespace aare

include/aare/ClusterFinderMT.hpp

@@ -0,0 +1,268 @@
+#pragma once
+#include <atomic>
+#include <cstdint>
+#include <memory>
+#include <thread>
+#include <vector>
+
+#include "aare/ClusterFinder.hpp"
+#include "aare/NDArray.hpp"
+#include "aare/ProducerConsumerQueue.hpp"
+
+namespace aare {
+
+enum class FrameType {
+    DATA,
+    PEDESTAL,
+};
+
+struct FrameWrapper {
+    FrameType type;
+    uint64_t frame_number;
+    NDArray<uint16_t, 2> data;
+};
+
+/**
+ * @brief ClusterFinderMT is a multi-threaded version of ClusterFinder. It uses
+ * a producer-consumer queue to distribute the frames to the threads. The
+ * clusters are collected in a single output queue.
+ * @tparam FRAME_TYPE type of the frame data
+ * @tparam PEDESTAL_TYPE type of the pedestal data
+ * @tparam CT type of the cluster data
+ */
+template <typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double,
+          typename CT = int32_t>
+class ClusterFinderMT {
+    size_t m_current_thread{0};
+    size_t m_n_threads{0};
+    using Finder = ClusterFinder<FRAME_TYPE, PEDESTAL_TYPE, CT>;
+    using InputQueue = ProducerConsumerQueue<FrameWrapper>;
+    using OutputQueue = ProducerConsumerQueue<ClusterVector<int>>;
+    std::vector<std::unique_ptr<InputQueue>> m_input_queues;
+    std::vector<std::unique_ptr<OutputQueue>> m_output_queues;
+
+    OutputQueue m_sink{1000}; // All clusters go into this queue
+
+    std::vector<std::unique_ptr<Finder>> m_cluster_finders;
+    std::vector<std::thread> m_threads;
+    std::thread m_collect_thread;
+    std::chrono::milliseconds m_default_wait{1};
+
+    std::atomic<bool> m_stop_requested{false};
+    std::atomic<bool> m_processing_threads_stopped{true};
+
+    /**
+     * @brief Function called by the processing threads. It reads the frames
+     * from the input queue and processes them.
+     */
+    void process(int thread_id) {
+        auto cf = m_cluster_finders[thread_id].get();
+        auto q = m_input_queues[thread_id].get();
+        bool realloc_same_capacity = true;
+
+        while (!m_stop_requested || !q->isEmpty()) {
+            if (FrameWrapper *frame = q->frontPtr(); frame != nullptr) {
+
+                switch (frame->type) {
+                case FrameType::DATA:
+                    cf->find_clusters(frame->data.view(), frame->frame_number);
+                    m_output_queues[thread_id]->write(cf->steal_clusters(realloc_same_capacity));
+                    break;
+
+                case FrameType::PEDESTAL:
+                    m_cluster_finders[thread_id]->push_pedestal_frame(
+                        frame->data.view());
+                    break;
+                }
+
+                // frame is processed now discard it
+                m_input_queues[thread_id]->popFront();
+            } else {
+                std::this_thread::sleep_for(m_default_wait);
+            }
+        }
+    }
+
+    /**
+     * @brief Collect all the clusters from the output queues and write them to
+     * the sink
+     */
+    void collect() {
+        bool empty = true;
+        while (!m_stop_requested || !empty || !m_processing_threads_stopped) {
+            empty = true;
+            for (auto &queue : m_output_queues) {
+                if (!queue->isEmpty()) {
+
+                    while (!m_sink.write(std::move(*queue->frontPtr()))) {
+                        std::this_thread::sleep_for(m_default_wait);
+                    }
+                    queue->popFront();
+                    empty = false;
+                }
+            }
+        }
+    }
+
+  public:
+    /**
+     * @brief Construct a new ClusterFinderMT object
+     * @param image_size size of the image
+     * @param cluster_size size of the cluster
+     * @param nSigma number of sigma above the pedestal to consider a photon
+     * @param capacity initial capacity of the cluster vector. Should match
+     * expected number of clusters in a frame per frame.
+     * @param n_threads number of threads to use
+     */
+    ClusterFinderMT(Shape<2> image_size, Shape<2> cluster_size,
+                    PEDESTAL_TYPE nSigma = 5.0, size_t capacity = 2000,
+                    size_t n_threads = 3)
+        : m_n_threads(n_threads) {
+        for (size_t i = 0; i < n_threads; i++) {
+            m_cluster_finders.push_back(
+                std::make_unique<ClusterFinder<FRAME_TYPE, PEDESTAL_TYPE, CT>>(
+                    image_size, cluster_size, nSigma, capacity));
+        }
+        for (size_t i = 0; i < n_threads; i++) {
+            m_input_queues.emplace_back(std::make_unique<InputQueue>(200));
+            m_output_queues.emplace_back(std::make_unique<OutputQueue>(200));
+        }
+        //TODO! Should we start automatically?
+        start();
+    }
+
+    /**
+     * @brief Return the sink queue where all the clusters are collected
+     * @warning You need to empty this queue otherwise the cluster finder will wait forever
+     */
+    ProducerConsumerQueue<ClusterVector<int>> *sink() { return &m_sink; }
+
+    /**
+     * @brief Start all processing threads
+     */
+    void start() {
+        m_processing_threads_stopped = false;
+        m_stop_requested = false;
+
+        for (size_t i = 0; i < m_n_threads; i++) {
+            m_threads.push_back(
+                std::thread(&ClusterFinderMT::process, this, i));
+        }
+
+        m_collect_thread = std::thread(&ClusterFinderMT::collect, this);
+    }
+
+    /**
+     * @brief Stop all processing threads
+     */
+    void stop() {
+        m_stop_requested = true;
+
+        for (auto &thread : m_threads) {
+            thread.join();
+        }
+        m_threads.clear();
+
+        m_processing_threads_stopped = true;
+        m_collect_thread.join();
+    }
+
+    /**
+     * @brief Wait for all the queues to be empty. Mostly used for timing tests.
+     */
+    void sync() {
+        for (auto &q : m_input_queues) {
+            while (!q->isEmpty()) {
+                std::this_thread::sleep_for(m_default_wait);
+            }
+        }
+        for (auto &q : m_output_queues) {
+            while (!q->isEmpty()) {
+                std::this_thread::sleep_for(m_default_wait);
+            }
+        }
+        while (!m_sink.isEmpty()) {
+            std::this_thread::sleep_for(m_default_wait);
+        }
+    }
+
+    /**
+     * @brief Push a pedestal frame to all the cluster finders. The frames is
+     * expected to be dark. No photon finding is done. Just pedestal update.
+     */
+    void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {
+        FrameWrapper fw{FrameType::PEDESTAL, 0,
+                        NDArray(frame)}; // TODO! copies the data!
+
+        for (auto &queue : m_input_queues) {
+            while (!queue->write(fw)) {
+                std::this_thread::sleep_for(m_default_wait);
+            }
+        }
+    }
+
+    /**
+     * @brief Push the frame to the queue of the next available thread. Function
+     * returns once the frame is in a queue.
+     * @note Spin locks with a default wait if the queue is full.
+     */
+    void find_clusters(NDView<FRAME_TYPE, 2> frame, uint64_t frame_number = 0) {
+        FrameWrapper fw{FrameType::DATA, frame_number,
+                        NDArray(frame)}; // TODO! copies the data!
+        while (!m_input_queues[m_current_thread % m_n_threads]->write(fw)) {
+            std::this_thread::sleep_for(m_default_wait);
+        }
+        m_current_thread++;
+    }
+
+    void clear_pedestal() {
+        if (!m_processing_threads_stopped) {
+            throw std::runtime_error("ClusterFinderMT is still running");
+        }
+        for (auto &cf : m_cluster_finders) {
+            cf->clear_pedestal();
+        }
+    }
+
+    /**
+     * @brief Return the pedestal currently used by the cluster finder
+     * @param thread_index index of the thread
+     */
+    auto pedestal(size_t thread_index = 0) {
+        if (m_cluster_finders.empty()) {
+            throw std::runtime_error("No cluster finders available");
+        }
+        if (!m_processing_threads_stopped) {
+            throw std::runtime_error("ClusterFinderMT is still running");
+        }
+        if (thread_index >= m_cluster_finders.size()) {
+            throw std::runtime_error("Thread index out of range");
+        }
+        return m_cluster_finders[thread_index]->pedestal();
+    }
+
+    /**
+     * @brief Return the noise currently used by the cluster finder
+     * @param thread_index index of the thread
+     */
+    auto noise(size_t thread_index = 0) {
+        if (m_cluster_finders.empty()) {
+            throw std::runtime_error("No cluster finders available");
+        }
+        if (!m_processing_threads_stopped) {
+            throw std::runtime_error("ClusterFinderMT is still running");
+        }
+        if (thread_index >= m_cluster_finders.size()) {
+            throw std::runtime_error("Thread index out of range");
+        }
+        return m_cluster_finders[thread_index]->noise();
+    }
+
+    // void push(FrameWrapper&& frame) {
+    //     //TODO! need to loop until we are successful
+    //     auto rc = m_input_queue.write(std::move(frame));
+    //     fmt::print("pushed frame {}\n", rc);
+    // }
+};
+
+} // namespace aare

include/aare/ClusterVector.hpp

@@ -0,0 +1,275 @@
+#pragma once
+#include <algorithm>
+#include <array>
+#include <cstddef>
+#include <cstdint>
+#include <numeric>
+#include <vector>
+
+#include <fmt/core.h>
+
+namespace aare {
+
+/**
+ * @brief ClusterVector is a container for clusters of various sizes. It uses a
+ * contiguous memory buffer to store the clusters. It is templated on the data
+ * type and the coordinate type of the clusters.
+ * @note push_back can invalidate pointers to elements in the container
+ * @warning ClusterVector is currently move only to catch unintended copies, but
+ * this might change since there are probably use cases where copying is needed.
+ * @tparam T data type of the pixels in the cluster
+ * @tparam CoordType data type of the x and y coordinates of the cluster
+ * (normally int16_t)
+ */
+template <typename T, typename CoordType = int16_t> class ClusterVector {
+    using value_type = T;
+    size_t m_cluster_size_x;
+    size_t m_cluster_size_y;
+    std::byte *m_data{};
+    size_t m_size{0};
+    size_t m_capacity;
+    uint64_t m_frame_number{0}; // TODO! Check frame number size and type
+    /*
+    Format string used in the python bindings to create a numpy
+    array from the buffer
+    = - native byte order
+    h - short
+    d - double
+    i - int
+    */
+    constexpr static char m_fmt_base[] = "=h:x:\nh:y:\n({},{}){}:data:";
+
+  public:
+    /**
+     * @brief Construct a new ClusterVector object
+     * @param cluster_size_x size of the cluster in x direction
+     * @param cluster_size_y size of the cluster in y direction
+     * @param capacity initial capacity of the buffer in number of clusters
+     * @param frame_number frame number of the clusters. Default is 0, which is
+     * also used to indicate that the clusters come from many frames
+     */
+    ClusterVector(size_t cluster_size_x = 3, size_t cluster_size_y = 3,
+                  size_t capacity = 1024, uint64_t frame_number = 0)
+        : m_cluster_size_x(cluster_size_x), m_cluster_size_y(cluster_size_y),
+          m_capacity(capacity), m_frame_number(frame_number) {
+        allocate_buffer(capacity);
+    }
+
+    ~ClusterVector() { delete[] m_data; }
+
+    // Move constructor
+    ClusterVector(ClusterVector &&other) noexcept
+        : m_cluster_size_x(other.m_cluster_size_x),
+          m_cluster_size_y(other.m_cluster_size_y), m_data(other.m_data),
+          m_size(other.m_size), m_capacity(other.m_capacity),
+          m_frame_number(other.m_frame_number) {
+        other.m_data = nullptr;
+        other.m_size = 0;
+        other.m_capacity = 0;
+    }
+
+    // Move assignment operator
+    ClusterVector &operator=(ClusterVector &&other) noexcept {
+        if (this != &other) {
+            delete[] m_data;
+            m_cluster_size_x = other.m_cluster_size_x;
+            m_cluster_size_y = other.m_cluster_size_y;
+            m_data = other.m_data;
+            m_size = other.m_size;
+            m_capacity = other.m_capacity;
+            m_frame_number = other.m_frame_number;
+            other.m_data = nullptr;
+            other.m_size = 0;
+            other.m_capacity = 0;
+            other.m_frame_number = 0;
+        }
+        return *this;
+    }
+
+    /**
+     * @brief Reserve space for at least capacity clusters
+     * @param capacity number of clusters to reserve space for
+     * @note If capacity is less than the current capacity, the function does
+     * nothing.
+     */
+    void reserve(size_t capacity) {
+        if (capacity > m_capacity) {
+            allocate_buffer(capacity);
+        }
+    }
+
+    /**
+     * @brief Add a cluster to the vector
+     * @param x x-coordinate of the cluster
+     * @param y y-coordinate of the cluster
+     * @param data pointer to the data of the cluster
+     * @warning The data pointer must point to a buffer of size cluster_size_x *
+     * cluster_size_y * sizeof(T)
+     */
+    void push_back(CoordType x, CoordType y, const std::byte *data) {
+        if (m_size == m_capacity) {
+            allocate_buffer(m_capacity * 2);
+        }
+        std::byte *ptr = element_ptr(m_size);
+        *reinterpret_cast<CoordType *>(ptr) = x;
+        ptr += sizeof(CoordType);
+        *reinterpret_cast<CoordType *>(ptr) = y;
+        ptr += sizeof(CoordType);
+
+        std::copy(data, data + m_cluster_size_x * m_cluster_size_y * sizeof(T),
+                  ptr);
+        m_size++;
+    }
+    ClusterVector &operator+=(const ClusterVector &other) {
+        if (m_size + other.m_size > m_capacity) {
+            allocate_buffer(m_capacity + other.m_size);
+        }
+        std::copy(other.m_data, other.m_data + other.m_size * item_size(),
+                  m_data + m_size * item_size());
+        m_size += other.m_size;
+        return *this;
+    }
+
+    /**
+     * @brief Sum the pixels in each cluster
+     * @return std::vector<T> vector of sums for each cluster
+     */
+    std::vector<T> sum() {
+        std::vector<T> sums(m_size);
+        const size_t stride = item_size();
+        const size_t n_pixels = m_cluster_size_x * m_cluster_size_y;
+        std::byte *ptr = m_data + 2 * sizeof(CoordType); // skip x and y
+
+        for (size_t i = 0; i < m_size; i++) {
+            sums[i] =
+                std::accumulate(reinterpret_cast<T *>(ptr),
+                                reinterpret_cast<T *>(ptr) + n_pixels, T{});
+            ptr += stride;
+        }
+        return sums;
+    }
+
+    /**
+     * @brief Return the maximum sum of the 2x2 subclusters in each cluster
+     * @return std::vector<T> vector of sums for each cluster
+     * @throws std::runtime_error if the cluster size is not 3x3
+     * @warning Only 3x3 clusters are supported for the 2x2 sum.
+     */
+    std::vector<T> sum_2x2() {
+        std::vector<T> sums(m_size);
+        const size_t stride = item_size();
+
+        if (m_cluster_size_x != 3 || m_cluster_size_y != 3) {
+            throw std::runtime_error(
+                "Only 3x3 clusters are supported for the 2x2 sum.");
+        }
+        std::byte *ptr = m_data + 2 * sizeof(CoordType); // skip x and y
+
+        for (size_t i = 0; i < m_size; i++) {
+            std::array<T, 4> total;
+            auto T_ptr = reinterpret_cast<T *>(ptr);
+            total[0] = T_ptr[0] + T_ptr[1] + T_ptr[3] + T_ptr[4];
+            total[1] = T_ptr[1] + T_ptr[2] + T_ptr[4] + T_ptr[5];
+            total[2] = T_ptr[3] + T_ptr[4] + T_ptr[6] + T_ptr[7];
+            total[3] = T_ptr[4] + T_ptr[5] + T_ptr[7] + T_ptr[8];
+
+            sums[i] = *std::max_element(total.begin(), total.end());
+            ptr += stride;
+        }
+
+        return sums;
+    }
+
+    /**
+     * @brief Return the number of clusters in the vector
+     */
+    size_t size() const { return m_size; }
+
+    /**
+     * @brief Return the capacity of the buffer in number of clusters. This is
+     * the number of clusters that can be stored in the current buffer without
+     * reallocation.
+     */
+    size_t capacity() const { return m_capacity; }
+
+    /**
+     * @brief Return the size in bytes of a single cluster
+     */
+    size_t item_size() const {
+        return 2 * sizeof(CoordType) +
+               m_cluster_size_x * m_cluster_size_y * sizeof(T);
+    }
+
+    /**
+     * @brief Return the offset in bytes for the i-th cluster
+     */
+    size_t element_offset(size_t i) const { return item_size() * i; }
+
+    /**
+     * @brief Return a pointer to the i-th cluster
+     */
+    std::byte *element_ptr(size_t i) { return m_data + element_offset(i); }
+
+    /**
+     * @brief Return a pointer to the i-th cluster
+     */
+    const std::byte *element_ptr(size_t i) const {
+        return m_data + element_offset(i);
+    }
+
+    size_t cluster_size_x() const { return m_cluster_size_x; }
+    size_t cluster_size_y() const { return m_cluster_size_y; }
+
+    std::byte *data() { return m_data; }
+    std::byte const *data() const { return m_data; }
+
+    /**
+     * @brief Return a reference to the i-th cluster casted to type V
+     * @tparam V type of the cluster
+     */
+    template <typename V> V &at(size_t i) {
+        return *reinterpret_cast<V *>(element_ptr(i));
+    }
+
+    const std::string_view fmt_base() const {
+        // TODO! how do we match on coord_t?
+        return m_fmt_base;
+    }
+
+    /**
+     * @brief Return the frame number of the clusters. 0 is used to indicate
+     * that the clusters come from many frames
+     */
+    uint64_t frame_number() const { return m_frame_number; }
+
+    void set_frame_number(uint64_t frame_number) {
+        m_frame_number = frame_number;
+    }
+
+    /**
+     * @brief Resize the vector to contain new_size clusters. If new_size is
+     * greater than the current capacity, a new buffer is allocated. If the size
+     * is smaller no memory is freed, size is just updated.
+     * @param new_size new size of the vector
+     * @warning The additional clusters are not initialized
+     */
+    void resize(size_t new_size) {
+        // TODO! Should we initialize the new clusters?
+        if (new_size > m_capacity) {
+            allocate_buffer(new_size);
+        }
+        m_size = new_size;
+    }
+
+  private:
+    void allocate_buffer(size_t new_capacity) {
+        size_t num_bytes = item_size() * new_capacity;
+        std::byte *new_data = new std::byte[num_bytes]{};
+        std::copy(m_data, m_data + item_size() * m_size, new_data);
+        delete[] m_data;
+        m_data = new_data;
+        m_capacity = new_capacity;
+    }
+};
+
+} // namespace aare

include/aare/CtbRawFile.hpp

@@ -0,0 +1,41 @@
+#pragma once
+
+#include "aare/FileInterface.hpp"
+#include "aare/RawMasterFile.hpp"
+#include "aare/Frame.hpp"
+
+#include <filesystem>
+#include <fstream>
+
+namespace aare{
+
+
+class CtbRawFile{
+    RawMasterFile m_master;
+    std::ifstream m_file;
+    size_t m_current_frame{0};
+    size_t m_current_subfile{0};
+    size_t m_num_subfiles{0};
+public:
+    CtbRawFile(const std::filesystem::path &fname);
+
+    void read_into(std::byte *image_buf, DetectorHeader* header = nullptr);
+    void seek(size_t frame_index);          //!< seek to the given frame index
+    size_t tell() const;                    //!< get the frame index of the file pointer
+
+    // in the specific class we can expose more functionality
+
+    size_t image_size_in_bytes() const;
+    size_t frames_in_file() const;
+
+    RawMasterFile master() const;
+private:
+    void find_subfiles();
+    size_t sub_file_index(size_t frame_index) const {
+        return frame_index / m_master.max_frames_per_file();
+    }
+    void open_data_file(size_t subfile_index);
+
+};
+
+}

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,66 @@
+#pragma once
+#include "aare/FileInterface.hpp"
+#include <memory>
+
+namespace aare {
+
+/**
+ * @brief RAII File class for reading, and in the future potentially writing 
+ * image files in various formats. Minimal generic interface. For specail fuctions
+ * plase use the RawFile or NumpyFile classes directly.
+ * Wraps FileInterface to abstract the underlying file format
+ * @note **frame_number** refers the the frame number sent by the detector while **frame_index**
+ * is the position of the frame in the file
+ */
+class File {
+    std::unique_ptr<FileInterface> file_impl;
+
+  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="r", const FileConfig &cfg = {});
+    
+     /**Since the object is responsible for managing the file we disable copy construction */
+    File(File const &other) = delete;     
+
+    /**The same goes for copy assignment */
+    File& operator=(File const &other) = delete;
+
+    File(File &&other) noexcept;
+    File& operator=(File &&other) noexcept;
+    ~File() = default;
+
+    // void close();                             //!< close the file
+      
+    Frame read_frame();                       //!< read one frame from the file at the current position
+    Frame read_frame(size_t frame_index);     //!< read one frame at the position given by frame number
+    std::vector<Frame> read_n(size_t n_frames); //!< read n_frames from the file at the current position
+
+    void read_into(std::byte *image_buf);
+    void read_into(std::byte *image_buf, size_t n_frames);
+    
+    size_t frame_number();                    //!< get the frame number at the current position
+    size_t frame_number(size_t frame_index);  //!< get the frame number at the given frame index
+    size_t bytes_per_frame() const;           
+    size_t pixels_per_frame() const;    
+    size_t bytes_per_pixel() const;  
+    size_t bitdepth() const;
+    void seek(size_t frame_index);          //!< seek to the given frame index
+    size_t tell() const;                     //!< get the frame index of the file pointer
+    size_t total_frames() const;
+    size_t rows() const;
+    size_t cols() const;
+
+    DetectorType detector_type() const;
+
+    
+};
+
+} // namespace aare

include/aare/FileInterface.hpp

@@ -0,0 +1,161 @@
+#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  one frame from the file at the current position
+     * @return Frame
+     */
+    virtual Frame read_frame() = 0;
+
+    /**
+     * @brief read one frame from the file at the given frame number
+     * @param frame_number frame number to read
+     * @return frame
+     */
+    virtual Frame read_frame(size_t frame_number) = 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_n(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;
+
+
+    virtual DetectorType detector_type() const = 0;
+
+    // function to query the data type of the file
+    /*virtual DataType dtype = 0; */
+
+    virtual ~FileInterface() = default;
+
+  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/Fit.hpp

@@ -0,0 +1,92 @@
+#pragma once
+
+#include <cmath>
+#include <fmt/core.h>
+#include <vector>
+
+#include "aare/NDArray.hpp"
+
+namespace aare {
+
+namespace func {
+double gaus(const double x, const double *par);
+NDArray<double, 1> gaus(NDView<double, 1> x, NDView<double, 1> par);
+
+double pol1(const double x, const double *par);
+NDArray<double, 1> pol1(NDView<double, 1> x, NDView<double, 1> par);
+
+} // namespace func
+
+
+/**
+ * @brief Estimate the initial parameters for a Gaussian fit
+ */
+std::array<double, 3> gaus_init_par(const NDView<double, 1> x, const NDView<double, 1> y);
+
+std::array<double, 2> pol1_init_par(const NDView<double, 1> x, const NDView<double, 1> y);
+
+static constexpr int DEFAULT_NUM_THREADS = 4;
+
+/**
+ * @brief Fit a 1D Gaussian to data.
+ * @param data data to fit
+ * @param x x values
+ */
+NDArray<double, 1> fit_gaus(NDView<double, 1> x, NDView<double, 1> y);
+
+
+/**
+ * @brief Fit a 1D Gaussian to each pixel. Data layout [row, col, values]
+ * @param x x values
+ * @param y y vales, layout [row, col, values]
+ * @param n_threads number of threads to use
+ */
+
+NDArray<double, 3> fit_gaus(NDView<double, 1> x, NDView<double, 3> y,
+                            int n_threads = DEFAULT_NUM_THREADS);
+
+
+
+
+/**
+ * @brief Fit a 1D Gaussian with error estimates
+ * @param x x values
+ * @param y y vales, layout [row, col, values]
+ * @param y_err error in y, layout [row, col, values]
+ * @param par_out output parameters
+ * @param par_err_out output error parameters
+ */
+void fit_gaus(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
+              NDView<double, 1> par_out, NDView<double, 1> par_err_out,
+              double& chi2);
+
+/**
+ * @brief Fit a 1D Gaussian to each pixel with error estimates. Data layout
+ * [row, col, values]
+ * @param x x values
+ * @param y y vales, layout [row, col, values]
+ * @param y_err error in y, layout [row, col, values]
+ * @param par_out output parameters, layout [row, col, values]
+ * @param par_err_out output parameter errors, layout [row, col, values]
+ * @param n_threads number of threads to use
+ */
+void fit_gaus(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
+              NDView<double, 3> par_out, NDView<double, 3> par_err_out, NDView<double, 2> chi2_out,
+              int n_threads = DEFAULT_NUM_THREADS
+              );
+
+NDArray<double, 1> fit_pol1(NDView<double, 1> x, NDView<double, 1> y);
+
+NDArray<double, 3> fit_pol1(NDView<double, 1> x, NDView<double, 3> y,
+                            int n_threads = DEFAULT_NUM_THREADS);
+
+void fit_pol1(NDView<double, 1> x, NDView<double, 1> y, NDView<double, 1> y_err,
+              NDView<double, 1> par_out, NDView<double, 1> par_err_out, double& chi2);
+
+// TODO! not sure we need to offer the different version in C++
+void fit_pol1(NDView<double, 1> x, NDView<double, 3> y, NDView<double, 3> y_err,
+              NDView<double, 3> par_out, NDView<double, 3> par_err_out,NDView<double, 2> chi2_out,
+              int n_threads = DEFAULT_NUM_THREADS);
+
+
+} // namespace aare

include/aare/Frame.hpp

@@ -0,0 +1,124 @@
+#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. Not much more than a
+ * pointer and some info. Limited interface to accept frames from many sources.
+ */
+class Frame {
+    uint32_t m_rows;
+    uint32_t m_cols;
+    Dtype m_dtype;
+    std::byte *m_data;
+    //TODO! Add frame number?
+
+  public:
+    /**
+     * @brief Construct a new Frame
+     * @param rows number of rows
+     * @param cols number of columns
+     * @param dtype data type of the pixels
+     * @note the data is initialized to zero
+     */
+    Frame(uint32_t rows, uint32_t cols, Dtype dtype);
+
+    /**
+     * @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 dtype data type of the pixels
+     */
+    Frame(const std::byte *bytes, uint32_t rows, uint32_t cols, Dtype dtype);
+    ~Frame(){ delete[] m_data; };
+
+    /** @warning Copy is disabled to ensure performance when passing
+     * frames around. Can discuss enabling it.
+     *
+     */
+    Frame &operator=(const Frame &other) = delete;
+    Frame(const Frame &other) = delete;
+
+    // enable move
+    Frame &operator=(Frame &&other) noexcept;
+    Frame(Frame &&other) noexcept;
+
+
+    Frame clone() const; //<- Explicit copy
+
+    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;
+
+    /**
+     * @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
+     * @warning The user should cast the pointer to the appropriate type. Think
+     * twice if this is the function you want to use.
+     */
+    std::byte *pixel_ptr(uint32_t row, uint32_t col) const;
+
+    /**
+     * @brief Set the pixel at the given row and column to the given value
+     * @tparam T type of the value
+     * @param row row index
+     * @param col column index
+     * @param data value to set
+     */
+    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(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");
+        }
+        //TODO! add tests then reimplement using pixel_ptr
+        T data;
+        std::memcpy(&data, m_data + (row * m_cols + col) * m_dtype.bytes(),
+                    m_dtype.bytes());
+        return data;
+    }
+    /**
+     * @brief Return an NDView of the frame. This is the preferred way to access
+     * data in the frame.
+     * 
+     * @tparam T type of the pixels
+     * @return NDView<T, 2> 
+     */
+    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);
+    }
+
+    /** 
+     * @brief Copy the frame data into a new NDArray. This is a deep copy.
+     */
+    template <typename T> NDArray<T> image() {
+        return NDArray<T>(this->view<T>());
+    }
+};
+
+} // namespace aare

include/aare/NDArray.hpp

@@ -0,0 +1,450 @@
+#pragma once
+/*
+Container holding image data, or a time series of image data in contigious
+memory.
+
+
+TODO! Add expression templates for operators
+
+*/
+#include "aare/ArrayExpr.hpp"
+#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 ArrayExpr<NDArray<T, Ndim>, Ndim> {
+    std::array<int64_t, Ndim> shape_;
+    std::array<int64_t, Ndim> strides_;
+    size_t size_{};
+    T *data_;
+
+  public:
+    /**
+     * @brief Default constructor. Will construct an empty NDArray.
+     *
+     */
+    NDArray() : shape_(), strides_(c_strides<Ndim>(shape_)), data_(nullptr) {};
+
+    /**
+     * @brief Construct a new NDArray object with a given shape.
+     * @note The data is uninitialized.
+     *
+     * @param shape shape of the new NDArray
+     */
+    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_]) {}
+
+
+    /**
+     * @brief Construct a new NDArray object with a shape and value.
+     *
+     * @param shape shape of the new array
+     * @param value value to initialize the array with
+     */
+    NDArray(std::array<int64_t, Ndim> shape, T value) : NDArray(shape) {
+        this->operator=(value);
+    }
+
+    /**
+     * @brief Construct a new NDArray object from a NDView.
+     * @note The data is copied from the view to the NDArray.
+     *
+     * @param v view of data to initialize the NDArray with
+     */
+    explicit NDArray(const NDView<T, Ndim> v) : NDArray(v.shape()) {
+        std::copy(v.begin(), v.end(), begin());
+    }
+
+    template<size_t Size>
+    NDArray(const std::array<T, Size>& arr) : NDArray<T,1>({Size}) {
+        std::copy(arr.begin(), arr.end(), begin());
+    }
+
+    // Move constructor
+    NDArray(NDArray &&other) noexcept
+        : shape_(other.shape_), strides_(c_strides<Ndim>(shape_)),
+          size_(other.size_), data_(other.data_) {
+        other.reset(); // TODO! is this necessary?
+
+    }
+
+    // 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_);
+    }
+
+    // Conversion operator from array expression to array
+    template <typename E>
+    NDArray(ArrayExpr<E, Ndim> &&expr) : NDArray(expr.shape()) {
+        for (size_t i = 0; i < size_; ++i) {
+            data_[i] = expr[i];
+        }
+    }
+
+    ~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);
+
+    //Write directly to the data array, or create a new one
+    template<size_t Size>
+    NDArray<T,1>& operator=(const std::array<T,Size> &other){
+        if(Size != size_){
+            delete[] data_;
+            size_ = Size;
+            data_ = new T[size_];
+        }
+        for (size_t i = 0; i < Size; ++i) {
+            data_[i] = other[i];
+        }
+        return *this;
+    }
+
+    // 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...)];
+    }
+
+    // TODO! is int the right type for index?
+    T &operator()(int64_t i) { return data_[i]; }
+    const T &operator()(int64_t i) const { return data_[i]; }
+
+    T &operator[](int64_t i) { return data_[i]; }
+    const T &operator[](int64_t i) const { return data_[i]; }
+
+    T *data() { return data_; }
+    std::byte *buffer() { return reinterpret_cast<std::byte *>(data_); }
+    size_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;
+    }
+
+    /**
+     * @brief Create a view of the NDArray.
+     *
+     * @return NDView<T, Ndim>
+     */
+    NDView<T, Ndim> view() 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);
+    }
+};
+
+// 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<T, Ndim> &other) {
+    // check shape
+    if (shape_ == other.shape_) {
+        for (size_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<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<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 T &mask) {
+    for (auto it = begin(); it != end(); ++it)
+        *it &= mask;
+    return *this;
+}
+
+template <typename T, int64_t Ndim>
+NDArray<bool, Ndim> NDArray<T, Ndim>::operator>(const NDArray &other) {
+    if (shape_ == other.shape_) {
+        NDArray<bool, Ndim> 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>
+std::ostream &operator<<(std::ostream &os, const NDArray<T, Ndim> &arr) {
+    for (auto row = 0; row < arr.shape(0); ++row) {
+        for (auto col = 0; col < arr.shape(1); ++col) {
+            os << std::setw(3);
+            os << arr(row, col) << " ";
+        }
+        os << "\n";
+    }
+    return os;
+}
+
+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,187 @@
+#pragma once
+#include "aare/defs.hpp"
+#include "aare/ArrayExpr.hpp"
+
+#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 ArrayExpr<NDView<T, Ndim>, Ndim> {
+  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...)];
+    }
+
+    size_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 const *begin() const { return buffer_; }
+    T const *end() const { 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>()); }
+
+
+    template<size_t Size>
+    NDView& operator=(const std::array<T, Size> &arr) {
+        if(size() != arr.size())
+            throw std::runtime_error(LOCATION + "Array and NDView size mismatch");
+        std::copy(arr.begin(), arr.end(), begin());
+        return *this;
+    }
+
+    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() const { 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";
+    }
+}
+
+
+template <typename T, int64_t Ndim>
+std::ostream& operator <<(std::ostream& os, const NDView<T, Ndim>& arr){
+    for (auto row = 0; row < arr.shape(0); ++row) {
+        for (auto col = 0; col < arr.shape(1); ++col) {
+            os << std::setw(3);
+            os << arr(row, col) << " ";
+        }
+        os << "\n";
+    }
+    return os;
+}
+
+
+} // namespace aare

include/aare/NumpyFile.hpp

@@ -0,0 +1,119 @@
+#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_frame() override { return get_frame(this->current_frame++); }
+    Frame read_frame(size_t frame_number) override { return get_frame(frame_number); }
+
+    std::vector<Frame> read_n(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(); }
+
+    DetectorType detector_type() const override { return DetectorType::Unknown; }
+
+    /**
+     * @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{};
+
+    size_t m_cols;
+    size_t m_rows;
+    size_t m_bitdepth;
+
+    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,209 @@
+#pragma once
+#include "aare/Frame.hpp"
+#include "aare/NDArray.hpp"
+#include "aare/NDView.hpp"
+#include <cstddef>
+
+namespace aare {
+
+/**
+ * @brief Calculate the pedestal of a series of frames. Can be used as
+ * standalone but mostly used in the ClusterFinder.
+ *
+ * @tparam SUM_TYPE type of the sum
+ */
+template <typename SUM_TYPE = double> class Pedestal {
+    uint32_t m_rows;
+    uint32_t m_cols;
+
+    uint32_t m_samples;
+    NDArray<uint32_t, 2> m_cur_samples;
+    
+    //TODO! in case of int needs to be changed to uint64_t
+    NDArray<SUM_TYPE, 2> m_sum;
+    NDArray<SUM_TYPE, 2> m_sum2;
+
+    //Cache mean since it is used over and over in the ClusterFinder
+    //This optimization is related to the access pattern of the ClusterFinder
+    //Relies on having more reads than pushes to the pedestal
+    NDArray<SUM_TYPE, 2> m_mean; 
+
+  public:
+    Pedestal(uint32_t rows, uint32_t cols, uint32_t n_samples = 1000)
+        : m_rows(rows), m_cols(cols), 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})),
+          m_mean(NDArray<SUM_TYPE, 2>({rows, cols})) {
+        assert(rows > 0 && cols > 0 && n_samples > 0);
+        m_sum = 0;
+        m_sum2 = 0;
+        m_mean = 0;
+    }
+    ~Pedestal() = default;
+
+    NDArray<SUM_TYPE, 2> mean() {
+        return m_mean;
+    }
+
+    SUM_TYPE mean(const uint32_t row, const uint32_t col) const {
+        return m_mean(row, col);
+    }
+
+    SUM_TYPE std(const uint32_t row, const uint32_t col) const {
+        return std::sqrt(variance(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);
+    }
+
+    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> std() {
+        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) =
+                std(i / m_cols, i % m_cols);
+        }
+
+        return standard_deviation_array;
+    }
+
+    
+
+    void clear() {
+        m_sum = 0;
+        m_sum2 = 0;
+        m_cur_samples = 0;
+        m_mean = 0;
+    }
+
+    
+
+    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;
+        m_mean(row, col) = 0;
+    }
+    
+
+
+    template <typename T> void push(NDView<T, 2> frame) {
+        assert(frame.size() == m_rows * m_cols);
+
+        // TODO! move away from m_rows, m_cols
+        if (frame.shape() != std::array<int64_t, 2>{m_rows, m_cols}) {
+            throw std::runtime_error(
+                "Frame shape does not match pedestal shape");
+        }
+
+        for (size_t row = 0; row < m_rows; row++) {
+            for (size_t col = 0; col < m_cols; col++) {
+                push<T>(row, col, frame(row, col));
+            }
+        }
+    }
+
+    /**
+     * Push but don't update the cached mean. Speeds up the process
+     * when initializing the pedestal.
+     * 
+     */
+    template <typename T> void push_no_update(NDView<T, 2> frame) {
+        assert(frame.size() == m_rows * m_cols);
+
+        // TODO! move away from m_rows, m_cols
+        if (frame.shape() != std::array<int64_t, 2>{m_rows, m_cols}) {
+            throw std::runtime_error(
+                "Frame shape does not match pedestal shape");
+        }
+
+        for (size_t row = 0; row < m_rows; row++) {
+            for (size_t col = 0; col < m_cols; col++) {
+                push_no_update<T>(row, col, frame(row, col));
+            }
+        }
+    }
+
+
+
+
+    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
+    uint32_t rows() const { return m_rows; }
+    uint32_t cols() const { return m_cols; }
+    uint32_t n_samples() const { return m_samples; }
+    NDArray<uint32_t, 2> cur_samples() const { return m_cur_samples; }
+    NDArray<SUM_TYPE, 2> get_sum() const { return m_sum; }
+    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>
+    void push(const uint32_t row, const uint32_t col, const T val_) {
+        SUM_TYPE val = static_cast<SUM_TYPE>(val_);
+        if (m_cur_samples(row, col) < m_samples) {
+            m_sum(row, col) += val;
+            m_sum2(row, col) += val * val;
+            m_cur_samples(row, col)++;
+        } else {
+            m_sum(row, col) += val - m_sum(row, col) / m_samples;
+            m_sum2(row, col) += val * val - m_sum2(row, col) / m_samples;
+        }
+        //Since we just did a push we know that m_cur_samples(row, col) is at least 1
+        m_mean(row, col) = m_sum(row, col) / m_cur_samples(row, col);
+    }
+
+    template <typename T>
+    void push_no_update(const uint32_t row, const uint32_t col, const T val_) {
+        SUM_TYPE val = static_cast<SUM_TYPE>(val_);
+        if (m_cur_samples(row, col) < m_samples) {
+            m_sum(row, col) += val;
+            m_sum2(row, col) += val * val;
+            m_cur_samples(row, col)++;
+        } else {
+            m_sum(row, col) += val - m_sum(row, col) / m_cur_samples(row, col);
+            m_sum2(row, col) += val * val - m_sum2(row, col) / m_cur_samples(row, col);
+        }
+    }
+
+    /**
+     * @brief Update the mean of the pedestal. This is used after having done
+     * push_no_update. It is not necessary to call this function after push.
+     */
+    void update_mean(){
+        m_mean = m_sum / m_cur_samples;
+    }
+
+    template<typename T>
+    void push_fast(const uint32_t row, const uint32_t col, const T val_){
+        //Assume we reached the steady state where all pixels have
+        //m_samples samples
+        SUM_TYPE val = static_cast<SUM_TYPE>(val_);
+        m_sum(row, col) += val - m_sum(row, col) / m_samples;
+        m_sum2(row, col) += val * val - m_sum2(row, col) / m_samples;
+        m_mean(row, col) = m_sum(row, col) / m_samples;
+    }
+
+};
+} // namespace aare

include/aare/PixelMap.hpp

@@ -0,0 +1,20 @@
+#pragma once
+
+#include "aare/defs.hpp"
+#include "aare/NDArray.hpp"
+
+namespace aare {
+
+NDArray<ssize_t, 2> GenerateMoench03PixelMap();
+NDArray<ssize_t, 2> GenerateMoench05PixelMap();
+NDArray<ssize_t, 2> GenerateMoench05PixelMap1g();
+NDArray<ssize_t, 2> GenerateMoench05PixelMapOld();
+
+//Matterhorn02
+NDArray<ssize_t, 2>GenerateMH02SingleCounterPixelMap();
+NDArray<ssize_t, 3> GenerateMH02FourCounterPixelMap();
+
+//Eiger 
+NDArray<ssize_t, 2>GenerateEigerFlipRowsPixelMap();
+
+} // namespace aare

include/aare/ProducerConsumerQueue.hpp

@@ -0,0 +1,203 @@
+/*
+ * Copyright (c) Meta Platforms, Inc. and affiliates.
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+// @author Bo Hu (bhu@fb.com)
+// @author Jordan DeLong (delong.j@fb.com)
+
+// Changes made by PSD Detector Group:
+// Copied: Line 34 constexpr std::size_t hardware_destructive_interference_size = 128; from folly/lang/Align.h
+// Changed extension to .hpp
+// Changed namespace to aare
+
+#pragma once
+
+#include <atomic>
+#include <cassert>
+#include <cstdlib>
+#include <memory>
+#include <stdexcept>
+#include <type_traits>
+#include <utility>
+
+constexpr std::size_t hardware_destructive_interference_size = 128;
+namespace aare {
+
+/*
+ * ProducerConsumerQueue is a one producer and one consumer queue
+ * without locks.
+ */
+template <class T> struct ProducerConsumerQueue {
+    typedef T value_type;
+
+    ProducerConsumerQueue(const ProducerConsumerQueue &) = delete;
+    ProducerConsumerQueue &operator=(const ProducerConsumerQueue &) = delete;
+
+    
+    ProducerConsumerQueue(ProducerConsumerQueue &&other){
+        size_ = other.size_;
+        records_ = other.records_;
+        other.records_ = nullptr;
+        readIndex_ = other.readIndex_.load(std::memory_order_acquire);
+        writeIndex_ = other.writeIndex_.load(std::memory_order_acquire);
+    }
+    ProducerConsumerQueue &operator=(ProducerConsumerQueue &&other){
+        size_ = other.size_;
+        records_ = other.records_;
+        other.records_ = nullptr;
+        readIndex_ = other.readIndex_.load(std::memory_order_acquire);
+        writeIndex_ = other.writeIndex_.load(std::memory_order_acquire);
+        return *this;
+    }
+    
+    
+    ProducerConsumerQueue():ProducerConsumerQueue(2){};
+    // size must be >= 2.
+    //
+    // Also, note that the number of usable slots in the queue at any
+    // given time is actually (size-1), so if you start with an empty queue,
+    // isFull() will return true after size-1 insertions.
+    explicit ProducerConsumerQueue(uint32_t size)
+        : size_(size), records_(static_cast<T *>(std::malloc(sizeof(T) * size))), readIndex_(0), writeIndex_(0) {
+        assert(size >= 2);
+        if (!records_) {
+            throw std::bad_alloc();
+        }
+    }
+
+    ~ProducerConsumerQueue() {
+        // We need to destruct anything that may still exist in our queue.
+        // (No real synchronization needed at destructor time: only one
+        // thread can be doing this.)
+        if (!std::is_trivially_destructible<T>::value) {
+            size_t readIndex = readIndex_;
+            size_t endIndex = writeIndex_;
+            while (readIndex != endIndex) {
+                records_[readIndex].~T();
+                if (++readIndex == size_) {
+                    readIndex = 0;
+                }
+            }
+        }
+
+        std::free(records_);
+    }
+
+    template <class... Args> bool write(Args &&...recordArgs) {
+        auto const currentWrite = writeIndex_.load(std::memory_order_relaxed);
+        auto nextRecord = currentWrite + 1;
+        if (nextRecord == size_) {
+            nextRecord = 0;
+        }
+        if (nextRecord != readIndex_.load(std::memory_order_acquire)) {
+            new (&records_[currentWrite]) T(std::forward<Args>(recordArgs)...);
+            writeIndex_.store(nextRecord, std::memory_order_release);
+            return true;
+        }
+
+        // queue is full
+        return false;
+    }
+
+    // move (or copy) the value at the front of the queue to given variable
+    bool read(T &record) {
+        auto const currentRead = readIndex_.load(std::memory_order_relaxed);
+        if (currentRead == writeIndex_.load(std::memory_order_acquire)) {
+            // queue is empty
+            return false;
+        }
+
+        auto nextRecord = currentRead + 1;
+        if (nextRecord == size_) {
+            nextRecord = 0;
+        }
+        record = std::move(records_[currentRead]);
+        records_[currentRead].~T();
+        readIndex_.store(nextRecord, std::memory_order_release);
+        return true;
+    }
+
+    // pointer to the value at the front of the queue (for use in-place) or
+    // nullptr if empty.
+    T *frontPtr() {
+        auto const currentRead = readIndex_.load(std::memory_order_relaxed);
+        if (currentRead == writeIndex_.load(std::memory_order_acquire)) {
+            // queue is empty
+            return nullptr;
+        }
+        return &records_[currentRead];
+    }
+
+    // queue must not be empty
+    void popFront() {
+        auto const currentRead = readIndex_.load(std::memory_order_relaxed);
+        assert(currentRead != writeIndex_.load(std::memory_order_acquire));
+
+        auto nextRecord = currentRead + 1;
+        if (nextRecord == size_) {
+            nextRecord = 0;
+        }
+        records_[currentRead].~T();
+        readIndex_.store(nextRecord, std::memory_order_release);
+    }
+
+    bool isEmpty() const {
+        return readIndex_.load(std::memory_order_acquire) == writeIndex_.load(std::memory_order_acquire);
+    }
+
+    bool isFull() const {
+        auto nextRecord = writeIndex_.load(std::memory_order_acquire) + 1;
+        if (nextRecord == size_) {
+            nextRecord = 0;
+        }
+        if (nextRecord != readIndex_.load(std::memory_order_acquire)) {
+            return false;
+        }
+        // queue is full
+        return true;
+    }
+
+    // * If called by consumer, then true size may be more (because producer may
+    //   be adding items concurrently).
+    // * If called by producer, then true size may be less (because consumer may
+    //   be removing items concurrently).
+    // * It is undefined to call this from any other thread.
+    size_t sizeGuess() const {
+        int ret = writeIndex_.load(std::memory_order_acquire) - readIndex_.load(std::memory_order_acquire);
+        if (ret < 0) {
+            ret += size_;
+        }
+        return ret;
+    }
+
+    // maximum number of items in the queue.
+    size_t capacity() const { return size_ - 1; }
+
+  private:
+    using AtomicIndex = std::atomic<unsigned int>;
+
+    char pad0_[hardware_destructive_interference_size];
+    // const uint32_t size_;
+    uint32_t size_;
+    // T *const records_;
+    T* records_;
+
+    alignas(hardware_destructive_interference_size) AtomicIndex readIndex_;
+    alignas(hardware_destructive_interference_size) AtomicIndex writeIndex_;
+
+    char pad1_[hardware_destructive_interference_size - sizeof(AtomicIndex)];
+};
+
+} // namespace aare

include/aare/RawFile.hpp

@@ -0,0 +1,126 @@
+#pragma once
+#include "aare/FileInterface.hpp"
+#include "aare/RawMasterFile.hpp"
+#include "aare/Frame.hpp"
+#include "aare/NDArray.hpp" //for pixel map
+#include "aare/RawSubFile.hpp"
+
+
+#include <optional>
+
+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 Class to read .raw files. The class will parse the master file
+ * to find the correct geometry for the frames.
+ * @note A more generic interface is available in the aare::File class.
+ * Consider using that unless you need raw file specific functionality.
+ */
+class RawFile : public FileInterface {
+    size_t n_subfiles{}; //f0,f1...fn
+    size_t n_subfile_parts{}; // d0,d1...dn
+    //TODO! move to vector of SubFile instead of pointers
+    std::vector<std::vector<RawSubFile *>> subfiles; //subfiles[f0,f1...fn][d0,d1...dn]
+    // std::vector<xy> positions;
+    
+    ModuleConfig cfg{0, 0};
+
+    RawMasterFile m_master;
+
+    size_t m_current_frame{};
+
+    // std::vector<ModuleGeometry> m_module_pixel_0;
+    // size_t m_rows{};
+    // size_t m_cols{};
+
+    DetectorGeometry m_geometry;
+
+  public:
+    /**
+     * @brief RawFile constructor
+     * @param fname path to the master file (.json)
+     * @param mode file mode (only "r" is supported at the moment)
+
+     */
+    RawFile(const std::filesystem::path &fname, const std::string &mode = "r");
+    virtual ~RawFile() override;
+
+    Frame read_frame() override;
+    Frame read_frame(size_t frame_number) override;
+    std::vector<Frame> read_n(size_t n_frames) override;
+    void read_into(std::byte *image_buf) override;
+    void read_into(std::byte *image_buf, size_t n_frames) override;
+
+    //TODO! do we need to adapt the API? 
+    void read_into(std::byte *image_buf, DetectorHeader *header);
+    void read_into(std::byte *image_buf, size_t n_frames, DetectorHeader *header);
+
+
+    size_t frame_number(size_t frame_index) override;
+    size_t bytes_per_frame() override;
+    size_t pixels_per_frame() override;
+    size_t bytes_per_pixel() const;
+    void seek(size_t frame_index) override;
+    size_t tell() override;
+    size_t total_frames() const override;
+    size_t rows() const override;
+    size_t cols() const override;
+    size_t bitdepth() const override;
+    xy geometry();
+    size_t n_mod() const;
+    
+    RawMasterFile master() const;
+
+
+
+    DetectorType detector_type() const override;
+
+  private:
+
+    /**
+     * @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, DetectorHeader *header = nullptr);
+
+
+    /**
+     * @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 read the header of the file
+     * @param fname path to the data subfile
+     * @return DetectorHeader
+     */
+    static DetectorHeader read_header(const std::filesystem::path &fname);
+
+    // void update_geometry_with_roi();
+    int find_number_of_subfiles();
+
+    void open_subfiles();
+    void find_geometry();
+};
+
+
+} // namespace aare

include/aare/RawMasterFile.hpp

@@ -0,0 +1,142 @@
+#pragma once
+#include "aare/defs.hpp"
+#include <filesystem>
+#include <fmt/format.h>
+#include <fstream>
+#include <optional>
+
+#include <nlohmann/json.hpp>
+using json = nlohmann::json;
+
+namespace aare {
+
+/**
+ * @brief Implementation used in RawMasterFile to parse the file name
+ */
+class RawFileNameComponents {
+    bool m_old_scheme{false};
+    std::filesystem::path m_base_path{};
+    std::string m_base_name{};
+    std::string m_ext{};
+    int m_file_index{}; // TODO! is this measurement_index?
+
+  public:
+    RawFileNameComponents(const std::filesystem::path &fname);
+
+    /// @brief Get the filename including path of the master file.
+    /// (i.e. what was passed in to the constructor))
+    std::filesystem::path master_fname() const;
+
+    /// @brief Get the filename including path of the data file.
+    /// @param mod_id module id run_d[module_id]_f0_0
+    /// @param file_id file id run_d0_f[file_id]_0
+    std::filesystem::path data_fname(size_t mod_id, size_t file_id) const;
+
+    const std::filesystem::path &base_path() const;
+    const std::string &base_name() const;
+    const std::string &ext() const;
+    int file_index() const;
+    void set_old_scheme(bool old_scheme);
+};
+
+class ScanParameters {
+    bool m_enabled = false;
+    std::string m_dac;
+    int m_start = 0;
+    int m_stop = 0;
+    int m_step = 0;
+    //TODO! add settleTime, requires string to time conversion
+
+  public:
+    ScanParameters(const std::string &par);
+    ScanParameters() = default;
+    ScanParameters(const ScanParameters &) = default;
+    ScanParameters &operator=(const ScanParameters &) = default;
+    ScanParameters(ScanParameters &&) = default;
+    int start() const;
+    int stop() const;
+    int step() const;
+    const std::string &dac() const;
+    bool enabled() const;
+    void increment_stop();
+};
+
+
+/**
+ * @brief Class for parsing a master file either in our .json format or the old
+ * .raw format
+ */
+class RawMasterFile {
+    RawFileNameComponents m_fnc;
+    std::string m_version;
+    DetectorType m_type;
+    TimingMode m_timing_mode;
+
+    size_t m_image_size_in_bytes{};
+    size_t m_frames_in_file{};
+    size_t m_total_frames_expected{};
+    size_t m_pixels_y{};
+    size_t m_pixels_x{};
+    size_t m_bitdepth{};
+
+    xy m_geometry{};
+
+    size_t m_max_frames_per_file{};
+    // uint32_t m_adc_mask{}; // TODO! implement reading
+    FrameDiscardPolicy m_frame_discard_policy{};
+    size_t m_frame_padding{};
+
+    // TODO! should these be bool?
+    uint8_t m_analog_flag{};
+    uint8_t m_digital_flag{};
+    uint8_t m_transceiver_flag{};
+
+    ScanParameters m_scan_parameters;
+
+    std::optional<size_t> m_analog_samples;
+    std::optional<size_t> m_digital_samples;
+    std::optional<size_t> m_transceiver_samples;
+    std::optional<size_t> m_number_of_rows;
+    std::optional<uint8_t> m_quad;
+
+    std::optional<ROI> m_roi;
+
+
+  public:
+    RawMasterFile(const std::filesystem::path &fpath);
+
+    std::filesystem::path data_fname(size_t mod_id, size_t file_id) const;
+
+    const std::string &version() const; //!< For example "7.2"
+    const DetectorType &detector_type() const;
+    const TimingMode &timing_mode() const;
+    size_t image_size_in_bytes() const;
+    size_t frames_in_file() const;
+    size_t pixels_y() const;
+    size_t pixels_x() const;
+    size_t max_frames_per_file() const;
+    size_t bitdepth() const;
+    size_t frame_padding() const;
+    const FrameDiscardPolicy &frame_discard_policy() const;
+
+    size_t total_frames_expected() const;
+    xy geometry() const;
+
+    std::optional<size_t> analog_samples() const;
+    std::optional<size_t> digital_samples() const;
+    std::optional<size_t> transceiver_samples() const;
+    std::optional<size_t> number_of_rows() const;
+    std::optional<uint8_t> quad() const;
+
+
+    std::optional<ROI> roi() const;
+
+
+    ScanParameters scan_parameters() const;
+
+  private:
+    void parse_json(const std::filesystem::path &fpath);
+    void parse_raw(const std::filesystem::path &fpath);
+};
+
+} // namespace aare

include/aare/RawSubFile.hpp

@@ -0,0 +1,75 @@
+#pragma once
+#include "aare/Frame.hpp"
+#include "aare/defs.hpp"
+
+#include <cstdint>
+#include <filesystem>
+#include <map>
+#include <optional>
+
+namespace aare {
+
+/**
+ * @brief Class to read a singe subfile written in .raw format. Used from RawFile to read
+ * the entire detector. Can be used directly to read part of the image.
+ */
+class RawSubFile {
+  protected:
+    std::ifstream m_file;
+    DetectorType m_detector_type;
+    size_t m_bitdepth;
+    std::filesystem::path m_fname;
+    size_t m_rows{};
+    size_t m_cols{};
+    size_t m_bytes_per_frame{};
+    size_t n_frames{};
+    uint32_t m_pos_row{};
+    uint32_t m_pos_col{};
+ 
+    
+    std::optional<NDArray<ssize_t, 2>> m_pixel_map;
+
+  public:
+    /**
+     * @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
+     */
+    RawSubFile(const std::filesystem::path &fname, DetectorType detector,
+               size_t rows, size_t cols, size_t bitdepth, uint32_t pos_row = 0, uint32_t pos_col = 0);
+
+    ~RawSubFile() = default;
+    /**
+     * @brief Seek to the given frame number
+     * @note Puts the file pointer at the start of the header, not the start of the data
+     * @param frame_index frame position in file to seek to
+     * @throws std::runtime_error if the frame number is out of range
+     */
+    void seek(size_t frame_index);
+    size_t tell();
+
+    void read_into(std::byte *image_buf, DetectorHeader *header = nullptr);
+    void get_part(std::byte *buffer, size_t frame_index);
+    
+    void read_header(DetectorHeader *header);
+    
+    size_t rows() const;
+    size_t cols() const;
+    
+    size_t frame_number(size_t frame_index);
+
+    size_t bytes_per_frame() const { return m_bytes_per_frame; }
+    size_t pixels_per_frame() const { return m_rows * m_cols; }
+    size_t bytes_per_pixel() const { return m_bitdepth / 8; }
+
+private:
+  template <typename T>
+  void read_with_map(std::byte *image_buf);
+
+};
+
+} // 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(Shape<2> 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 (size_t 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 (size_t i = 0; i != labeled_.size(); ++i) {
+        auto cl = labeled_(i);
+        if (cl != 0) {
+            auto it = child.find(cl);
+            while (it != child.end()) {
+                cl = it->second;
+                it = child.find(cl);
+                // do this once before doing the second pass?
+                // all values point to the final one...
+            }
+            labeled_(i) = cl;
+        }
+    }
+}
+
+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_map;
+    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_map[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_map)
+        hits.push_back(h.second);
+}
+
+} // namespace aare

include/aare/decode.hpp

@@ -0,0 +1,13 @@
+#pragma once
+
+#include <cstdint>
+#include <aare/NDView.hpp>
+namespace aare {
+
+
+uint16_t adc_sar_05_decode64to16(uint64_t input);
+uint16_t adc_sar_04_decode64to16(uint64_t input);
+void adc_sar_05_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> output);
+void adc_sar_04_decode64to16(NDView<uint64_t, 2> input, NDView<uint16_t,2> output);
+
+} // namespace aare

include/aare/defs.hpp

@@ -0,0 +1,262 @@
+#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__) + ":"
+
+
+
+#ifdef AARE_CUSTOM_ASSERT
+#define AARE_ASSERT(expr)\
+    if (expr)\
+        {}\
+    else\
+        aare::assert_failed(LOCATION + " Assertion failed: " + #expr + "\n");
+#else
+#define AARE_ASSERT(cond)\
+    do { (void)sizeof(cond); } while(0)
+#endif
+
+
+namespace aare {
+
+void assert_failed(const std::string &msg);
+
+
+class DynamicCluster {
+  public:
+    int cluster_sizeX;
+    int cluster_sizeY;
+    int16_t x;
+    int16_t y;
+    Dtype dt; // 4 bytes
+
+  private:
+    std::byte *m_data;
+
+  public:
+    DynamicCluster(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()]{};
+    }
+    DynamicCluster() : DynamicCluster(3, 3) {}
+    DynamicCluster(const DynamicCluster &other)
+        : DynamicCluster(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());
+    }
+    DynamicCluster &operator=(const DynamicCluster &other) {
+        if (this == &other)
+            return *this;
+        this->~DynamicCluster();
+        new (this) DynamicCluster(other);
+        return *this;
+    }
+    DynamicCluster(DynamicCluster &&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);
+    }
+    ~DynamicCluster() { 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, dt.bytes());
+    }
+
+    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 DetectorHeader {
+    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>;
+
+
+/**
+ * @brief Class to hold the geometry of a module. Where pixel 0 is located and the size of the module
+ */
+struct ModuleGeometry{
+    int origin_x{};
+    int origin_y{};
+    int height{};
+    int width{};
+    int row_index{};
+    int col_index{}; 
+};
+
+/**
+ * @brief Class to hold the geometry of a detector. Number of modules, their size and where pixel 0 
+ * for each module is located
+ */
+struct DetectorGeometry{
+    int modules_x{};
+    int modules_y{};
+    int pixels_x{};
+    int pixels_y{};
+    int module_gap_row{};
+    int module_gap_col{};
+    std::vector<ModuleGeometry> module_pixel_0;
+};
+
+struct ROI{
+    int64_t xmin{};
+    int64_t xmax{};
+    int64_t ymin{};
+    int64_t ymax{};
+  
+    int64_t height() const { return ymax - ymin; }
+    int64_t width() const { return xmax - xmin; }
+  };
+
+
+using dynamic_shape = std::vector<int64_t>;
+
+//TODO! Can we uniform enums between the libraries?
+
+/**
+ * @brief Enum class to identify different detectors. 
+ * The values are the same as in slsDetectorPackage
+ * Different spelling to avoid confusion with the slsDetectorPackage
+ */
+enum class DetectorType {
+    //Standard detectors match the enum values from slsDetectorPackage
+    Generic,
+    Eiger,
+    Gotthard,
+    Jungfrau,
+    ChipTestBoard,
+    Moench,
+    Mythen3,
+    Gotthard2,
+    Xilinx_ChipTestBoard,
+
+    //Additional detectors used for defining processing. Variants of the standard ones.
+    Moench03=100,
+    Moench03_old,
+    Unknown
+};
+
+enum class TimingMode { Auto, Trigger };
+enum class FrameDiscardPolicy { NoDiscard, Discard, DiscardPartial };
+
+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*/);
+
+template <> FrameDiscardPolicy StringTo(const std::string & /*mode*/);
+
+using DataTypeVariants = std::variant<uint16_t, uint32_t>;
+
+} // namespace aare

include/aare/geo_helpers.hpp

@@ -0,0 +1,16 @@
+#pragma once
+#include "aare/defs.hpp"
+#include "aare/RawMasterFile.hpp" //ROI refactor away
+namespace aare{
+
+/**
+ * @brief Update the detector geometry given a region of interest
+ * 
+ * @param geo 
+ * @param roi 
+ * @return DetectorGeometry 
+ */
+DetectorGeometry update_geometry_with_roi(DetectorGeometry geo, ROI roi);
+
+
+} // namespace aare

include/aare/utils/par.hpp

@@ -0,0 +1,18 @@
+#include <thread>
+#include <vector>
+#include <utility>
+
+namespace aare {
+
+    template<typename F>
+    void RunInParallel(F func, const std::vector<std::pair<int, int>>& tasks) {
+        // auto tasks = split_task(0, y.shape(0), n_threads);
+        std::vector<std::thread> threads;
+        for (auto &task : tasks) {
+            threads.push_back(std::thread(func, task.first, task.second));
+        }
+        for (auto &thread : threads) {
+            thread.join();
+        }
+    }
+} // namespace aare

include/aare/utils/task.hpp

@@ -0,0 +1,8 @@
+
+#include <utility>
+#include <vector>
+
+namespace aare {
+std::vector<std::pair<int, int>> split_task(int first, int last, int n_threads);
+
+} // namespace aare

patches/lmfit.patch

@@ -0,0 +1,13 @@
+diff --git a/lib/CMakeLists.txt b/lib/CMakeLists.txt
+index 4efb7ed..6533660 100644
+--- a/lib/CMakeLists.txt
++++ b/lib/CMakeLists.txt
+@@ -11,7 +11,7 @@ target_compile_definitions(${lib} PRIVATE "LMFIT_EXPORT") # for Windows DLL expo
+ 
+ target_include_directories(${lib}
+     PUBLIC
+-    $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/>
++    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/>
+     $<INSTALL_INTERFACE:include/>
+     )
+ 

pyproject.toml

@@ -0,0 +1,16 @@
+[build-system]
+requires = ["scikit-build-core>=0.10", "pybind11", "numpy"]
+build-backend = "scikit_build_core.build"
+
+[project]
+name = "aare"
+version = "2025.2.18"
+
+
+[tool.scikit-build]
+cmake.verbose = true
+
+[tool.scikit-build.cmake.define]
+AARE_PYTHON_BINDINGS = "ON"
+AARE_SYSTEM_LIBRARIES = "ON"
+AARE_INSTALL_PYTHONEXT = "ON"

python/CMakeLists.txt

@@ -0,0 +1,67 @@
+
+find_package (Python 3.10 COMPONENTS Interpreter Development REQUIRED)
+
+# Download or find pybind11 depending on configuration
+if(AARE_FETCH_PYBIND11)
+    FetchContent_Declare(
+        pybind11
+        GIT_REPOSITORY https://github.com/pybind/pybind11
+        GIT_TAG        v2.13.0
+    )
+    FetchContent_MakeAvailable(pybind11)
+else()
+    find_package(pybind11 2.13 REQUIRED)
+endif()
+
+# Add the compiled python extension
+pybind11_add_module(
+    _aare           # name of the module
+    src/module.cpp  # source file
+    )
+
+set_target_properties(_aare PROPERTIES
+    LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}
+)
+target_link_libraries(_aare PRIVATE aare_core aare_compiler_flags)
+
+# List of python files to be copied to the build directory
+set( PYTHON_FILES
+    aare/__init__.py
+    aare/CtbRawFile.py
+    aare/func.py
+    aare/RawFile.py
+    aare/transform.py
+    aare/ScanParameters.py
+    aare/utils.py
+)
+
+# Copy the python files to the build directory
+foreach(FILE ${PYTHON_FILES})
+    configure_file(${FILE} ${CMAKE_BINARY_DIR}/${FILE}  )
+endforeach(FILE ${PYTHON_FILES})
+
+set_target_properties(_aare PROPERTIES
+    LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/aare
+)
+
+set(PYTHON_EXAMPLES
+    examples/play.py
+    examples/fits.py
+)
+
+
+# Copy the python examples to the build directory
+foreach(FILE ${PYTHON_EXAMPLES})
+    configure_file(${FILE} ${CMAKE_BINARY_DIR}/${FILE}  )
+    message(STATUS "Copying ${FILE} to ${CMAKE_BINARY_DIR}/${FILE}")
+endforeach(FILE ${PYTHON_EXAMPLES})
+
+
+if(AARE_INSTALL_PYTHONEXT)
+    install(TARGETS _aare
+        EXPORT "${TARGETS_EXPORT_NAME}"
+        LIBRARY DESTINATION aare
+    )
+
+    install(FILES ${PYTHON_FILES} DESTINATION aare)
+endif()

python/aare/CtbRawFile.py

@@ -0,0 +1,191 @@
+
+from . import _aare
+import numpy as np
+from .ScanParameters import ScanParameters
+
+class CtbRawFile(_aare.CtbRawFile):
+    """File reader for the CTB raw file format.
+    
+    Args:  
+            fname (pathlib.Path | str): Path to the file to be read.
+            chunk_size (int): Number of frames to read at a time. Default is 1.
+            transform (function): Function to apply to the data after reading it. 
+                The function should take a numpy array of type uint8 and return one
+                or several numpy arrays.
+    """
+    def __init__(self, fname, chunk_size = 1, transform = None):
+        super().__init__(fname)
+        self._chunk_size = chunk_size
+        self._transform = transform
+
+
+    def read_frame(self, frame_index: int | None = None ) -> tuple:
+        """Read one frame from the file and then advance the file pointer.
+
+        .. note::
+
+            Uses the position of the file pointer :py:meth:`~CtbRawFile.tell` to determine
+            which frame to read unless frame_index is specified.
+
+        Args:
+            frame_index (int): If not None, seek to this frame before reading.
+
+        Returns:
+            tuple: header, data
+
+        Raises:
+            RuntimeError: If the file is at the end.
+        """
+        if frame_index is not None:
+            self.seek(frame_index)
+
+
+        header, data = super().read_frame()
+        if header.shape == (1,):
+            header = header[0]
+
+
+        if self._transform:
+            res = self._transform(data)
+            if isinstance(res, tuple):
+                return header, *res
+            else:
+                return header, res
+        else:
+            return header, data
+            
+    def read_n(self, n_frames:int) -> tuple:
+        """Read several frames from the file.
+
+        .. note::
+
+            Uses the position of the file pointer :py:meth:`~CtbRawFile.tell` to determine
+            where to start reading from.
+
+        If the number of frames requested is larger than the number of frames left in the file,
+        the function will read the remaining frames. If no frames are left in the file
+        a RuntimeError is raised.
+
+        Args:
+            n_frames (int): Number of frames to read.
+
+        Returns:
+            tuple: header, data
+
+        Raises:
+            RuntimeError: If EOF is reached.
+        """
+        # Calculate the number of frames to actually read
+        n_frames = min(n_frames, self.frames_in_file - self.tell())
+        if n_frames == 0:
+            raise RuntimeError("No frames left in file.")
+
+
+        # Do the first read to figure out what we have
+        tmp_header, tmp_data = self.read_frame()
+        
+        # Allocate arrays for
+        header = np.zeros(n_frames, dtype = tmp_header.dtype)
+        data = np.zeros((n_frames, *tmp_data.shape), dtype = tmp_data.dtype)
+        
+        # Copy the first frame
+        header[0] = tmp_header
+        data[0] = tmp_data
+
+        # Do the rest of the reading
+        for i in range(1, n_frames):
+            header[i], data[i] = self.read_frame()
+
+        return header, data
+
+    def read(self) -> tuple:
+        """Read the entire file.
+        Seeks to the beginning of the file before reading.
+
+        Returns:
+            tuple: header, data
+        """
+        self.seek(0)
+        return self.read_n(self.frames_in_file)
+
+    def seek(self, frame_index:int) -> None:
+        """Seek to a specific frame in the file.
+
+        Args:
+            frame_index (int): Frame position in file to seek to.
+        """
+        super().seek(frame_index)
+
+    def tell(self) -> int:
+        """Return the current frame position in the file.
+
+        Returns:
+            int: Frame position in file.
+        """
+        return super().tell()
+    
+    @property
+    def scan_parameters(self):
+        """Return the scan parameters.
+
+        Returns:
+            ScanParameters: Scan parameters.
+        """
+        return ScanParameters(self.master.scan_parameters)
+
+    @property
+    def master(self):
+        """Return the master file.
+
+        Returns:
+            RawMasterFile: Master file.
+        """
+        return super().master()
+
+    @property
+    def image_size_in_bytes(self) -> int:
+        """Return the size of the image in bytes.
+
+        Returns:
+            int: Size of image in bytes.
+        """
+        return super().image_size_in_bytes
+
+    def __len__(self) -> int:
+        """Return the number of frames in the file.
+
+        Returns:
+            int: Number of frames in file.
+        """
+        return super().frames_in_file
+    
+    @property
+    def frames_in_file(self) -> int:
+        """Return the number of frames in the file.
+
+        Returns:
+            int: Number of frames in file.
+        """
+        return super().frames_in_file
+
+    def __enter__(self):
+        return self
+    
+    def __exit__(self, exc_type, exc_value, traceback):
+        pass
+
+    def __iter__(self):
+        return self
+    
+    def __next__(self):
+        try:
+            if self._chunk_size == 1:
+                return self.read_frame()
+            else:
+                return self.read_n(self._chunk_size)
+            
+                
+        except RuntimeError:
+            # TODO! find a good way to check that we actually have the right exception
+            raise StopIteration
+

python/aare/RawFile.py

@@ -0,0 +1,66 @@
+from . import _aare
+import numpy as np
+from .ScanParameters import ScanParameters
+
+class RawFile(_aare.RawFile):
+    def __init__(self, fname, chunk_size = 1):
+        super().__init__(fname)
+        self._chunk_size = chunk_size
+
+
+    def read(self) -> tuple:
+        """Read the entire file.
+        Seeks to the beginning of the file before reading.
+
+        Returns:
+            tuple: header, data
+        """
+        self.seek(0)
+        return self.read_n(self.total_frames)
+
+    @property
+    def scan_parameters(self):
+        """Return the scan parameters.
+
+        Returns:
+            ScanParameters: Scan parameters.
+        """
+        return ScanParameters(self.master.scan_parameters)
+    
+    @property
+    def master(self):
+        """Return the master file.
+
+        Returns:
+            RawMasterFile: Master file.
+        """
+        return super().master
+    
+    def __len__(self) -> int:
+        """Return the number of frames in the file.
+
+        Returns:
+            int: Number of frames in file.
+        """
+        return super().frames_in_file
+
+    def __enter__(self):
+        return self
+    
+    def __exit__(self, exc_type, exc_value, traceback):
+        pass
+
+    def __iter__(self):
+        return self
+    
+    def __next__(self):
+        try:
+            if self._chunk_size == 1:
+                return self.read_frame()
+            else:
+                return self.read_n(self._chunk_size)
+            
+                
+        except RuntimeError:
+            # TODO! find a good way to check that we actually have the right exception
+            raise StopIteration

python/aare/ScanParameters.py

@@ -0,0 +1,16 @@
+from . import _aare
+
+class ScanParameters(_aare.ScanParameters):
+    def __init__(self, s):
+        super().__init__(s)
+        
+    def __iter__(self):
+        return [getattr(self, a) for a in ['start', 'stop', 'step']].__iter__()
+    
+    def __str__(self):
+        return f'ScanParameters({self.dac}: {self.start}, {self.stop}, {self.step})'
+    
+    def __repr__(self):
+        return self.__str__()
+    
+    

python/aare/__init__.py

@@ -0,0 +1,23 @@
+# Make the compiled classes that live in _aare available from aare.
+from . import _aare
+
+
+from ._aare import File, RawMasterFile, RawSubFile
+from ._aare import Pedestal_d, Pedestal_f, ClusterFinder, VarClusterFinder
+from ._aare import DetectorType
+from ._aare import ClusterFile
+from ._aare import hitmap
+
+from ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i
+
+from ._aare import fit_gaus, fit_pol1
+
+from .CtbRawFile import CtbRawFile
+from .RawFile import RawFile
+from .ScanParameters import ScanParameters
+
+from .utils import random_pixels, random_pixel, flat_list
+
+
+#make functions available in the top level API
+from .func import *

python/aare/func.py

@@ -0,0 +1 @@
+from ._aare import gaus, pol1

python/aare/transform.py

@@ -0,0 +1,58 @@
+import numpy as np
+from . import _aare
+
+
+class AdcSar04Transform64to16:
+    def __call__(self, data):
+        return _aare.adc_sar_04_decode64to16(data)
+
+class AdcSar05Transform64to16:
+    def __call__(self, data):
+        return _aare.adc_sar_05_decode64to16(data)
+
+class Moench05Transform:
+    #Could be moved to C++ without changing the interface
+    def __init__(self):
+        self.pixel_map = _aare.GenerateMoench05PixelMap()
+
+    def __call__(self, data):
+        return np.take(data.view(np.uint16), self.pixel_map)
+    
+
+class Moench05Transform1g:
+    #Could be moved to C++ without changing the interface
+    def __init__(self):
+        self.pixel_map = _aare.GenerateMoench05PixelMap1g()
+
+    def __call__(self, data):
+        return np.take(data.view(np.uint16), self.pixel_map)
+    
+
+class Moench05TransformOld:
+    #Could be moved to C++ without changing the interface
+    def __init__(self):
+        self.pixel_map = _aare.GenerateMoench05PixelMapOld()
+
+    def __call__(self, data):
+        return np.take(data.view(np.uint16), self.pixel_map)
+
+
+class Matterhorn02Transform:
+    def __init__(self):
+        self.pixel_map = _aare.GenerateMH02FourCounterPixelMap()
+
+    def __call__(self, data):
+        counters = int(data.size / 48**2 / 2)
+        if counters == 1:
+            return np.take(data.view(np.uint16), self.pixel_map[0])
+        else:
+            return np.take(data.view(np.uint16), self.pixel_map[0:counters])
+
+
+#on import generate the pixel maps to avoid doing it every time
+moench05 = Moench05Transform()
+moench05_1g = Moench05Transform1g()
+moench05_old = Moench05TransformOld()
+matterhorn02 = Matterhorn02Transform()
+adc_sar_04_64to16 = AdcSar04Transform64to16()
+adc_sar_05_64to16 = AdcSar05Transform64to16()

python/aare/utils.py

@@ -0,0 +1,27 @@
+import numpy as np
+
+def random_pixels(n_pixels, xmin=0, xmax=512, ymin=0, ymax=1024):
+    """Return a list of random pixels.
+    
+    Args:
+        n_pixels (int): Number of pixels to return.
+        rows (int): Number of rows in the image.
+        cols (int): Number of columns in the image.
+        
+    Returns:
+        list: List of (row, col) tuples.
+    """
+    return [(np.random.randint(ymin, ymax), np.random.randint(xmin, xmax)) for _ in range(n_pixels)]
+
+
+def random_pixel(xmin=0, xmax=512, ymin=0, ymax=1024):
+    """Return a random pixel.
+    
+    Returns:
+        tuple: (row, col)
+    """
+    return random_pixels(1, xmin, xmax, ymin, ymax)[0]
+
+def flat_list(xss):
+    """Flatten a list of lists."""
+    return [x for xs in xss for x in xs]

python/examples/fits.py

@@ -0,0 +1,79 @@
+import matplotlib.pyplot as plt
+import numpy as np
+from aare import fit_gaus, fit_pol1
+from aare import gaus, pol1
+
+textpm = f"±"  #
+textmu = f"μ"  #
+textsigma = f"σ"  #
+
+
+
+# ================================= Gauss fit =================================
+# Parameters
+mu = np.random.uniform(1, 100)  # Mean of Gaussian
+sigma = np.random.uniform(4, 20)  # Standard deviation
+num_points = 10000  # Number of points for smooth distribution
+noise_sigma = 100
+
+# Generate Gaussian distribution
+data = np.random.normal(mu, sigma, num_points)
+
+# Generate errors for each point
+errors = np.abs(np.random.normal(0, sigma, num_points))  # Errors with mean 0, std 0.5
+
+# Create subplot
+fig0, ax0 = plt.subplots(1, 1, num=0, figsize=(12, 8))
+
+x = np.histogram(data, bins=30)[1][:-1] + 0.05
+y = np.histogram(data, bins=30)[0]
+yerr = errors[:30]
+
+
+# Add the errors as error bars in the step plot
+ax0.errorbar(x, y, yerr=yerr, fmt=". ", capsize=5)
+ax0.grid()
+
+par, err = fit_gaus(x, y, yerr)
+print(par, err)
+
+x = np.linspace(x[0], x[-1], 1000)
+ax0.plot(x, gaus(x, par), marker="")
+ax0.set(xlabel="x", ylabel="Counts", title=f"A0 = {par[0]:0.2f}{textpm}{err[0]:0.2f}\n"
+                                           f"{textmu} = {par[1]:0.2f}{textpm}{err[1]:0.2f}\n"
+                                           f"{textsigma} = {par[2]:0.2f}{textpm}{err[2]:0.2f}\n"
+                                           f"(init: {textmu}: {mu:0.2f}, {textsigma}: {sigma:0.2f})")
+fig0.tight_layout()
+
+
+
+# ================================= pol1 fit =================================
+# Parameters
+n_points = 40
+
+# Generate random slope and intercept (origin)
+slope = np.random.uniform(-10, 10)  # Random slope between 0.5 and 2.0
+intercept = np.random.uniform(-10, 10)  # Random intercept between -10 and 10
+
+# Generate random x values
+x_values = np.random.uniform(-10, 10, n_points)
+
+# Calculate y values based on the linear function y = mx + b + error
+errors = np.abs(np.random.normal(0, np.random.uniform(1, 5), n_points))
+var_points = np.random.normal(0, np.random.uniform(0.1, 2), n_points)
+y_values = slope * x_values + intercept + var_points
+
+fig1, ax1 = plt.subplots(1, 1, num=1, figsize=(12, 8))
+ax1.errorbar(x_values, y_values, yerr=errors, fmt=". ", capsize=5)
+par, err = fit_pol1(x_values, y_values, errors)
+
+
+x = np.linspace(np.min(x_values), np.max(x_values), 1000)
+ax1.plot(x, pol1(x, par), marker="")
+ax1.set(xlabel="x", ylabel="y", title=f"a = {par[0]:0.2f}{textpm}{err[0]:0.2f}\n"
+                                      f"b = {par[1]:0.2f}{textpm}{err[1]:0.2f}\n"
+                                      f"(init: {slope:0.2f}, {intercept:0.2f})")
+fig1.tight_layout()
+
+plt.show()
+

python/examples/play.py

@@ -0,0 +1,50 @@
+import sys
+sys.path.append('/home/l_msdetect/erik/aare/build')
+
+#Our normal python imports
+from pathlib import Path
+import matplotlib.pyplot as plt
+import numpy as np
+import boost_histogram as bh
+import time
+
+
+import aare
+
+data = np.random.normal(10, 1, 1000)
+
+hist = bh.Histogram(bh.axis.Regular(10, 0, 20))
+hist.fill(data)
+
+
+x = hist.axes[0].centers
+y = hist.values()
+y_err = np.sqrt(y)+1
+res = aare.fit_gaus(x, y, y_err, chi2 = True)
+
+
+        
+t_elapsed = time.perf_counter()-t0
+print(f'Histogram filling took: {t_elapsed:.3f}s {total_clusters/t_elapsed/1e6:.3f}M clusters/s')
+
+histogram_data = hist3d.counts()
+x = hist3d.axes[2].edges[:-1]
+
+y = histogram_data[100,100,:]
+xx = np.linspace(x[0], x[-1])
+# fig, ax = plt.subplots()
+# ax.step(x, y, where = 'post')
+
+y_err = np.sqrt(y)
+y_err = np.zeros(y.size)
+y_err += 1
+
+# par = fit_gaus2(y,x, y_err)
+# ax.plot(xx, gaus(xx,par))
+# print(par)
+
+res = fit_gaus(y,x)
+res2 = fit_gaus(y,x, y_err)
+print(res)
+print(res2)
+

python/src/cluster.hpp

@@ -0,0 +1,201 @@
+#include "aare/ClusterCollector.hpp"
+#include "aare/ClusterFileSink.hpp"
+#include "aare/ClusterFinder.hpp"
+#include "aare/ClusterFinderMT.hpp"
+#include "aare/ClusterVector.hpp"
+#include "aare/NDView.hpp"
+#include "aare/Pedestal.hpp"
+#include "np_helper.hpp"
+
+#include <cstdint>
+#include <filesystem>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl_bind.h>
+
+namespace py = pybind11;
+using pd_type = double;
+
+template <typename T>
+void define_cluster_vector(py::module &m, const std::string &typestr) {
+    auto class_name = fmt::format("ClusterVector_{}", typestr);
+    py::class_<ClusterVector<T>>(m, class_name.c_str(), py::buffer_protocol())
+        .def(py::init<int, int>())
+        .def_property_readonly("size", &ClusterVector<T>::size)
+        .def("item_size", &ClusterVector<T>::item_size)
+        .def_property_readonly("fmt",
+                               [typestr](ClusterVector<T> &self) {
+                                   return fmt::format(
+                                       self.fmt_base(), self.cluster_size_x(),
+                                       self.cluster_size_y(), typestr);
+                               })
+        .def("sum",
+             [](ClusterVector<T> &self) {
+                 auto *vec = new std::vector<T>(self.sum());
+                 return return_vector(vec);
+             })
+        .def("sum_2x2", [](ClusterVector<T> &self) {
+            auto *vec = new std::vector<T>(self.sum_2x2());
+            return return_vector(vec);
+        })
+        .def_property_readonly("capacity", &ClusterVector<T>::capacity)
+        .def_property("frame_number", &ClusterVector<T>::frame_number,
+                      &ClusterVector<T>::set_frame_number)
+        .def_buffer([typestr](ClusterVector<T> &self) -> py::buffer_info {
+            return py::buffer_info(
+                self.data(),      /* Pointer to buffer */
+                self.item_size(), /* Size of one scalar */
+                fmt::format(self.fmt_base(), self.cluster_size_x(),
+                            self.cluster_size_y(),
+                            typestr), /* Format descriptor */
+                1,                    /* Number of dimensions */
+                {self.size()},        /* Buffer dimensions */
+                {self.item_size()}    /* Strides (in bytes) for each index */
+            );
+        });
+}
+
+void define_cluster_finder_mt_bindings(py::module &m) {
+    py::class_<ClusterFinderMT<uint16_t, pd_type>>(m, "ClusterFinderMT")
+        .def(py::init<Shape<2>, Shape<2>, pd_type, size_t, size_t>(),
+             py::arg("image_size"), py::arg("cluster_size"),
+             py::arg("n_sigma") = 5.0, py::arg("capacity") = 2048,
+             py::arg("n_threads") = 3)
+        .def("push_pedestal_frame",
+             [](ClusterFinderMT<uint16_t, pd_type> &self,
+                py::array_t<uint16_t> frame) {
+                 auto view = make_view_2d(frame);
+                 self.push_pedestal_frame(view);
+             })
+        .def(
+            "find_clusters",
+            [](ClusterFinderMT<uint16_t, pd_type> &self,
+               py::array_t<uint16_t> frame, uint64_t frame_number) {
+                auto view = make_view_2d(frame);
+                self.find_clusters(view, frame_number);
+                return;
+            },
+            py::arg(), py::arg("frame_number") = 0)
+        .def("clear_pedestal", &ClusterFinderMT<uint16_t, pd_type>::clear_pedestal)
+        .def("sync", &ClusterFinderMT<uint16_t, pd_type>::sync)
+        .def("stop", &ClusterFinderMT<uint16_t, pd_type>::stop)
+        .def("start", &ClusterFinderMT<uint16_t, pd_type>::start)
+        .def("pedestal",
+             [](ClusterFinderMT<uint16_t, pd_type> &self, size_t thread_index) {
+                 auto pd = new NDArray<pd_type, 2>{};
+                 *pd = self.pedestal(thread_index);
+                 return return_image_data(pd);
+             },py::arg("thread_index") = 0)
+        .def("noise",
+             [](ClusterFinderMT<uint16_t, pd_type> &self, size_t thread_index) {
+                 auto arr = new NDArray<pd_type, 2>{};
+                 *arr = self.noise(thread_index);
+                 return return_image_data(arr);
+             },py::arg("thread_index") = 0);
+}
+
+void define_cluster_collector_bindings(py::module &m) {
+    py::class_<ClusterCollector>(m, "ClusterCollector")
+        .def(py::init<ClusterFinderMT<uint16_t, double, int32_t> *>())
+        .def("stop", &ClusterCollector::stop)
+        .def(
+            "steal_clusters",
+            [](ClusterCollector &self) {
+                auto v =
+                    new std::vector<ClusterVector<int>>(self.steal_clusters());
+                return v;
+            },
+            py::return_value_policy::take_ownership);
+}
+
+void define_cluster_file_sink_bindings(py::module &m) {
+    py::class_<ClusterFileSink>(m, "ClusterFileSink")
+        .def(py::init<ClusterFinderMT<uint16_t, double, int32_t> *,
+                      const std::filesystem::path &>())
+        .def("stop", &ClusterFileSink::stop);
+}
+
+void define_cluster_finder_bindings(py::module &m) {
+    py::class_<ClusterFinder<uint16_t, pd_type>>(m, "ClusterFinder")
+        .def(py::init<Shape<2>, Shape<2>, pd_type, size_t>(),
+             py::arg("image_size"), py::arg("cluster_size"),
+             py::arg("n_sigma") = 5.0, py::arg("capacity") = 1'000'000)
+        .def("push_pedestal_frame",
+             [](ClusterFinder<uint16_t, pd_type> &self,
+                py::array_t<uint16_t> frame) {
+                 auto view = make_view_2d(frame);
+                 self.push_pedestal_frame(view);
+             })
+        .def("clear_pedestal", &ClusterFinder<uint16_t, pd_type>::clear_pedestal)
+        .def_property_readonly("pedestal",
+                               [](ClusterFinder<uint16_t, pd_type> &self) {
+                                   auto pd = new NDArray<pd_type, 2>{};
+                                   *pd = self.pedestal();
+                                   return return_image_data(pd);
+                               })
+        .def_property_readonly("noise",
+                               [](ClusterFinder<uint16_t, pd_type> &self) {
+                                   auto arr = new NDArray<pd_type, 2>{};
+                                   *arr = self.noise();
+                                   return return_image_data(arr);
+                               })
+        .def(
+            "steal_clusters",
+            [](ClusterFinder<uint16_t, pd_type> &self,
+               bool realloc_same_capacity) {
+                auto v = new ClusterVector<int>(
+                    self.steal_clusters(realloc_same_capacity));
+                return v;
+            },
+            py::arg("realloc_same_capacity") = false)
+        .def(
+            "find_clusters",
+            [](ClusterFinder<uint16_t, pd_type> &self,
+               py::array_t<uint16_t> frame, uint64_t frame_number) {
+                auto view = make_view_2d(frame);
+                self.find_clusters(view, frame_number);
+                return;
+            },
+            py::arg(), py::arg("frame_number") = 0);
+
+    m.def("hitmap",
+          [](std::array<size_t, 2> image_size, ClusterVector<int32_t> &cv) {
+              py::array_t<int32_t> hitmap(image_size);
+              auto r = hitmap.mutable_unchecked<2>();
+
+              // Initialize hitmap to 0
+              for (py::ssize_t i = 0; i < r.shape(0); i++)
+                  for (py::ssize_t j = 0; j < r.shape(1); j++)
+                      r(i, j) = 0;
+
+              size_t stride = cv.item_size();
+              auto ptr = cv.data();
+              for (size_t i = 0; i < cv.size(); i++) {
+                  auto x = *reinterpret_cast<int16_t *>(ptr);
+                  auto y = *reinterpret_cast<int16_t *>(ptr + sizeof(int16_t));
+                  r(y, x) += 1;
+                  ptr += stride;
+              }
+              return hitmap;
+          });
+    define_cluster_vector<int>(m, "i");
+    define_cluster_vector<double>(m, "d");
+    define_cluster_vector<float>(m, "f");
+
+    py::class_<DynamicCluster>(m, "DynamicCluster", py::buffer_protocol())
+        .def(py::init<int, int, Dtype>())
+        .def("size", &DynamicCluster::size)
+        .def("begin", &DynamicCluster::begin)
+        .def("end", &DynamicCluster::end)
+        .def_readwrite("x", &DynamicCluster::x)
+        .def_readwrite("y", &DynamicCluster::y)
+        .def_buffer([](DynamicCluster &c) -> py::buffer_info {
+            return py::buffer_info(c.data(), c.dt.bytes(), c.dt.format_descr(),
+                                   1, {c.size()}, {c.dt.bytes()});
+        })
+
+        .def("__repr__", [](const DynamicCluster &a) {
+            return "<DynamicCluster: x: " + std::to_string(a.x) +
+                   ", y: " + std::to_string(a.y) + ">";
+        });
+}

python/src/cluster_file.hpp

@@ -0,0 +1,72 @@
+#include "aare/ClusterFile.hpp"
+#include "aare/defs.hpp"
+
+#include <cstdint>
+#include <filesystem>
+#include <pybind11/iostream.h>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl/filesystem.h>
+#include <string>
+
+//Disable warnings for unused parameters, as we ignore some
+//in the __exit__ method
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wunused-parameter"
+
+
+namespace py = pybind11;
+using namespace ::aare;
+
+void define_cluster_file_io_bindings(py::module &m) {
+    PYBIND11_NUMPY_DTYPE(Cluster3x3, x, y, data);
+
+    py::class_<ClusterFile>(m, "ClusterFile")
+        .def(py::init<const std::filesystem::path &, size_t,
+                      const std::string &>(),
+             py::arg(), py::arg("chunk_size") = 1000, py::arg("mode") = "r")
+        .def("read_clusters",
+             [](ClusterFile &self, size_t n_clusters) {
+                auto v = new ClusterVector<int32_t>(self.read_clusters(n_clusters));
+                return v;
+             },py::return_value_policy::take_ownership)
+        .def("read_frame",
+             [](ClusterFile &self) {
+                auto v = new ClusterVector<int32_t>(self.read_frame());
+                return v;
+             })
+        .def("write_frame", &ClusterFile::write_frame)
+        // .def("read_cluster_with_cut",
+        //      [](ClusterFile &self, size_t n_clusters,
+        //         py::array_t<double> noise_map, int nx, int ny) {
+        //          auto view = make_view_2d(noise_map);
+        //          auto *vec =
+        //              new std::vector<Cluster3x3>(self.read_cluster_with_cut(
+        //                  n_clusters, view.data(), nx, ny));
+        //          return return_vector(vec);
+        //      })
+        .def("__enter__", [](ClusterFile &self) { return &self; })
+        .def("__exit__",
+             [](ClusterFile &self,
+                const std::optional<pybind11::type> &exc_type,
+                const std::optional<pybind11::object> &exc_value,
+                const std::optional<pybind11::object> &traceback) {
+                 self.close();
+             })
+        .def("__iter__", [](ClusterFile &self) { return &self; })
+        .def("__next__", [](ClusterFile &self) {
+            auto v = new ClusterVector<int32_t>(self.read_clusters(self.chunk_size()));
+            if (v->size() == 0) {
+                throw py::stop_iteration();
+            }
+            return v;
+        });
+
+        m.def("calculate_eta2", []( aare::ClusterVector<int32_t> &clusters) {
+            auto eta2 = new NDArray<double, 2>(calculate_eta2(clusters));
+            return return_image_data(eta2);
+        });
+}
+
+#pragma GCC diagnostic pop

python/src/ctb_raw_file.hpp

@@ -0,0 +1,99 @@
+
+#include "aare/CtbRawFile.hpp"
+#include "aare/File.hpp"
+#include "aare/Frame.hpp"
+#include "aare/RawFile.hpp"
+#include "aare/RawMasterFile.hpp"
+#include "aare/RawSubFile.hpp"
+
+#include "aare/defs.hpp"
+#include "aare/decode.hpp"
+// #include "aare/fClusterFileV2.hpp"
+
+#include <cstdint>
+#include <filesystem>
+#include <pybind11/iostream.h>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl/filesystem.h>
+#include <string>
+
+namespace py = pybind11;
+using namespace ::aare;
+
+void define_ctb_raw_file_io_bindings(py::module &m) {
+
+m.def("adc_sar_05_decode64to16", [](py::array_t<uint8_t> input) {
+
+    
+    if(input.ndim() != 2){
+        throw std::runtime_error("Only 2D arrays are supported at this moment"); 
+    }
+
+    //Create a 2D output array with the same shape as the input
+    std::vector<ssize_t> shape{input.shape(0), input.shape(1)/8};
+    py::array_t<uint16_t> output(shape);
+
+    //Create a view of the input and output arrays
+    NDView<uint64_t, 2> input_view(reinterpret_cast<uint64_t*>(input.mutable_data()), {output.shape(0), output.shape(1)});
+    NDView<uint16_t, 2> output_view(output.mutable_data(), {output.shape(0), output.shape(1)});
+
+    adc_sar_05_decode64to16(input_view, output_view);
+
+    return output;
+});
+
+
+m.def("adc_sar_04_decode64to16", [](py::array_t<uint8_t> input) {
+
+    
+    if(input.ndim() != 2){
+        throw std::runtime_error("Only 2D arrays are supported at this moment"); 
+    }
+
+    //Create a 2D output array with the same shape as the input
+    std::vector<ssize_t> shape{input.shape(0), input.shape(1)/8};
+    py::array_t<uint16_t> output(shape);
+
+    //Create a view of the input and output arrays
+    NDView<uint64_t, 2> input_view(reinterpret_cast<uint64_t*>(input.mutable_data()), {output.shape(0), output.shape(1)});
+    NDView<uint16_t, 2> output_view(output.mutable_data(), {output.shape(0), output.shape(1)});
+
+    adc_sar_04_decode64to16(input_view, output_view);
+
+    return output;
+});
+
+ py::class_<CtbRawFile>(m, "CtbRawFile")
+        .def(py::init<const std::filesystem::path &>())
+        .def("read_frame",
+             [](CtbRawFile &self) {
+                 size_t image_size = self.image_size_in_bytes();
+                 py::array image;
+                 std::vector<ssize_t> shape;
+                 shape.reserve(2);
+                 shape.push_back(1);
+                 shape.push_back(image_size);
+
+                 py::array_t<DetectorHeader> header(1);
+
+                 // always read bytes
+                 image = py::array_t<uint8_t>(shape);
+
+                 self.read_into(
+                     reinterpret_cast<std::byte *>(image.mutable_data()),
+                     header.mutable_data());
+
+                 return py::make_tuple(header, image);
+             })
+        .def("seek", &CtbRawFile::seek)
+        .def("tell", &CtbRawFile::tell)
+        .def("master", &CtbRawFile::master)
+
+        .def_property_readonly("image_size_in_bytes",
+                               &CtbRawFile::image_size_in_bytes)
+
+        .def_property_readonly("frames_in_file", &CtbRawFile::frames_in_file);
+
+}

python/src/file.hpp

@@ -0,0 +1,295 @@
+#include "aare/CtbRawFile.hpp"
+#include "aare/File.hpp"
+#include "aare/Frame.hpp"
+#include "aare/RawFile.hpp"
+#include "aare/RawMasterFile.hpp"
+#include "aare/RawSubFile.hpp"
+
+#include "aare/defs.hpp"
+// #include "aare/fClusterFileV2.hpp"
+
+#include <cstdint>
+#include <filesystem>
+#include <pybind11/iostream.h>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl/filesystem.h>
+#include <string>
+
+namespace py = pybind11;
+using namespace ::aare;
+
+//Disable warnings for unused parameters, as we ignore some
+//in the __exit__ method
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wunused-parameter"
+
+void define_file_io_bindings(py::module &m) {
+
+
+    py::enum_<DetectorType>(m, "DetectorType")
+        .value("Jungfrau", DetectorType::Jungfrau)
+        .value("Eiger", DetectorType::Eiger)
+        .value("Mythen3", DetectorType::Mythen3)
+        .value("Moench", DetectorType::Moench)
+        .value("Moench03", DetectorType::Moench03)
+        .value("Moench03_old", DetectorType::Moench03_old)
+        .value("ChipTestBoard", DetectorType::ChipTestBoard)
+        .value("Unknown", DetectorType::Unknown);
+
+
+    PYBIND11_NUMPY_DTYPE(DetectorHeader, frameNumber, expLength, packetNumber,
+                         bunchId, timestamp, modId, row, column, reserved,
+                         debug, roundRNumber, detType, version, packetMask);
+
+   
+
+    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("frame_number", py::overload_cast<>(&File::frame_number))
+        .def("frame_number", py::overload_cast<size_t>(&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("bytes_per_pixel", &File::bytes_per_pixel)
+        .def_property_readonly(
+            "detector_type",
+            [](File &self) { return ToString(self.detector_type()); })
+        .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;
+             })
+        .def("read_frame",
+             [](File &self, size_t frame_number) {
+                 self.seek(frame_number);
+                 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;
+             })
+        .def("read_n", [](File &self, size_t n_frames) {
+            //adjust for actual frames left in the file
+            n_frames = std::min(n_frames, self.total_frames()-self.tell());
+            if(n_frames == 0){
+                throw std::runtime_error("No frames left in file");
+            }
+            std::vector<size_t> shape{n_frames, self.rows(), self.cols()};
+            
+            py::array image;
+            const uint8_t item_size = self.bytes_per_pixel();
+            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()),
+                           n_frames);
+            return image;
+        })
+        .def("__enter__", [](File &self) { return &self; })
+        .def("__exit__",
+             [](File &self,
+                const std::optional<pybind11::type> &exc_type,
+                const std::optional<pybind11::object> &exc_value,
+                const std::optional<pybind11::object> &traceback) {
+                //  self.close();
+             })
+        .def("__iter__", [](File &self) { return &self; })
+        .def("__next__", [](File &self) {
+
+            try{
+                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;
+            }catch(std::runtime_error &e){
+                throw py::stop_iteration();
+            }
+        });
+
+
+    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_<ScanParameters>(m, "ScanParameters")
+        .def(py::init<const std::string &>())
+        .def(py::init<const ScanParameters &>())
+
+        .def_property_readonly("enabled", &ScanParameters::enabled)
+        .def_property_readonly("dac", &ScanParameters::dac)
+        .def_property_readonly("start", &ScanParameters::start)
+        .def_property_readonly("stop", &ScanParameters::stop)
+        .def_property_readonly("step", &ScanParameters::step);
+
+
+    py::class_<ROI>(m, "ROI")
+        .def(py::init<>())
+        .def_readwrite("xmin", &ROI::xmin)
+        .def_readwrite("xmax", &ROI::xmax)
+        .def_readwrite("ymin", &ROI::ymin)
+        .def_readwrite("ymax", &ROI::ymax)
+        .def("__str__", [](const ROI& self){
+            return fmt::format("ROI: xmin: {} xmax: {} ymin: {} ymax: {}", self.xmin, self.xmax, self.ymin, self.ymax);
+        })
+        .def("__repr__", [](const ROI& self){
+            return fmt::format("<ROI: xmin: {} xmax: {} ymin: {} ymax: {}>", self.xmin, self.xmax, self.ymin, self.ymax);
+        })
+        .def("__iter__", [](const ROI &self) {
+            return py::make_iterator(&self.xmin, &self.ymax+1); //NOLINT
+        });
+
+    
+
+
+    py::class_<RawSubFile>(m, "RawSubFile")
+        .def(py::init<const std::filesystem::path &, DetectorType, size_t,
+                      size_t, size_t>())
+        .def_property_readonly("bytes_per_frame", &RawSubFile::bytes_per_frame)
+        .def_property_readonly("pixels_per_frame",
+                               &RawSubFile::pixels_per_frame)
+        .def("seek", &RawSubFile::seek)
+        .def("tell", &RawSubFile::tell)
+        .def_property_readonly("rows", &RawSubFile::rows)
+        .def_property_readonly("cols", &RawSubFile::cols)
+        .def("read_frame",
+             [](RawSubFile &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);
+                 }
+                 fmt::print("item_size: {} rows: {} cols: {}\n", item_size, self.rows(), self.cols());
+                 self.read_into(
+                     reinterpret_cast<std::byte *>(image.mutable_data()));
+                 return image;
+             });
+
+#pragma GCC diagnostic pop
+    // 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<DynamicCluster> &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/fit.hpp

@@ -0,0 +1,250 @@
+#include <cstdint>
+#include <filesystem>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl_bind.h>
+
+#include "aare/Fit.hpp"
+
+namespace py = pybind11;
+using namespace pybind11::literals;
+
+
+void define_fit_bindings(py::module &m) {
+
+    // TODO! Evaluate without converting to double
+    m.def(
+        "gaus",
+        [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
+           py::array_t<double, py::array::c_style | py::array::forcecast> par) {
+            auto x_view = make_view_1d(x);
+            auto par_view = make_view_1d(par);
+            auto y = new NDArray<double, 1>{aare::func::gaus(x_view, par_view)};
+            return return_image_data(y);
+        },
+        R"(
+        Evaluate a 1D Gaussian function for all points in x using parameters par.
+
+        Parameters
+        ----------
+        x : array_like
+            The points at which to evaluate the Gaussian function.
+        par : array_like
+            The parameters of the Gaussian function. The first element is the amplitude, the second element is the mean, and the third element is the standard deviation.
+        )",
+        py::arg("x"), py::arg("par"));
+
+    m.def(
+        "pol1",
+        [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
+           py::array_t<double, py::array::c_style | py::array::forcecast> par) {
+            auto x_view = make_view_1d(x);
+            auto par_view = make_view_1d(par);
+            auto y = new NDArray<double, 1>{aare::func::pol1(x_view, par_view)};
+            return return_image_data(y);
+        },
+        R"(
+        Evaluate a 1D polynomial function for all points in x using parameters par. (p0+p1*x)
+        
+        Parameters
+        ----------
+        x : array_like
+            The points at which to evaluate the polynomial function.
+        par : array_like
+            The parameters of the polynomial function. The first element is the intercept, and the second element is the slope.    
+        )",
+        py::arg("x"), py::arg("par"));
+
+
+    m.def(
+        "fit_gaus",
+        [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
+           py::array_t<double, py::array::c_style | py::array::forcecast> y,
+           int n_threads) {
+            if (y.ndim() == 3) {
+                auto par = new NDArray<double, 3>{};
+                auto y_view = make_view_3d(y);
+                auto x_view = make_view_1d(x);
+                *par = aare::fit_gaus(x_view, y_view, n_threads);
+                return return_image_data(par);
+            } else if (y.ndim() == 1) {
+                auto par = new NDArray<double, 1>{};
+                auto y_view = make_view_1d(y);
+                auto x_view = make_view_1d(x);
+                *par = aare::fit_gaus(x_view, y_view);
+                return return_image_data(par);
+            } else {
+                throw std::runtime_error("Data must be 1D or 3D");
+            }
+        },
+        R"(
+
+Fit a 1D Gaussian to data.
+
+Parameters
+----------
+x : array_like
+    The x values.
+y : array_like
+    The y values.
+n_threads : int, optional
+    The number of threads to use. Default is 4.
+)",
+        py::arg("x"), py::arg("y"), py::arg("n_threads") = 4);
+
+    m.def(
+        "fit_gaus",
+        [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
+           py::array_t<double, py::array::c_style | py::array::forcecast> y,
+           py::array_t<double, py::array::c_style | py::array::forcecast> y_err,
+           int n_threads) {
+
+            if (y.ndim() == 3) {
+                // Allocate memory for the output
+                // Need to have pointers to allow python to manage
+                // the memory
+                auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 3});
+                auto par_err =
+                    new NDArray<double, 3>({y.shape(0), y.shape(1), 3});
+                auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
+
+                // Make views of the numpy arrays
+                auto y_view = make_view_3d(y);
+                auto y_view_err = make_view_3d(y_err);
+                auto x_view = make_view_1d(x);
+
+                aare::fit_gaus(x_view, y_view, y_view_err, par->view(),
+                               par_err->view(), chi2->view(), n_threads);
+
+                return py::dict("par"_a = return_image_data(par),
+                                "par_err"_a = return_image_data(par_err),
+                                "chi2"_a = return_image_data(chi2),
+                                "Ndf"_a = y.shape(2) - 3);
+            } else if (y.ndim() == 1) {
+                // Allocate memory for the output
+                // Need to have pointers to allow python to manage
+                // the memory
+                auto par = new NDArray<double, 1>({3});
+                auto par_err = new NDArray<double, 1>({3});
+
+                // Decode the numpy arrays
+                auto y_view = make_view_1d(y);
+                auto y_view_err = make_view_1d(y_err);
+                auto x_view = make_view_1d(x);
+
+
+                double chi2 = 0;
+                aare::fit_gaus(x_view, y_view, y_view_err, par->view(),
+                               par_err->view(), chi2);
+
+                return py::dict("par"_a = return_image_data(par),
+                                "par_err"_a = return_image_data(par_err),
+                                "chi2"_a = chi2, "Ndf"_a = y.size() - 3);
+
+            } else {
+                throw std::runtime_error("Data must be 1D or 3D");
+            }
+        },
+        R"(
+
+Fit a 1D Gaussian to data with error estimates.
+
+Parameters
+----------
+x : array_like
+    The x values.
+y : array_like
+    The y values.
+y_err : array_like
+    The error in the y values.
+n_threads : int, optional
+    The number of threads to use. Default is 4.
+)",
+        py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
+
+    m.def(
+        "fit_pol1",
+        [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
+           py::array_t<double, py::array::c_style | py::array::forcecast> y,
+           int n_threads) {
+            if (y.ndim() == 3) {
+                auto par = new NDArray<double, 3>{};
+
+                auto x_view = make_view_1d(x);
+                auto y_view = make_view_3d(y);
+                *par = aare::fit_pol1(x_view, y_view, n_threads);
+                return return_image_data(par);
+            } else if (y.ndim() == 1) {
+                auto par = new NDArray<double, 1>{};
+                auto x_view = make_view_1d(x);
+                auto y_view = make_view_1d(y);
+                *par = aare::fit_pol1(x_view, y_view);
+                return return_image_data(par);
+            } else {
+                throw std::runtime_error("Data must be 1D or 3D");
+            }
+        },
+        py::arg("x"), py::arg("y"), py::arg("n_threads") = 4);
+
+    m.def(
+        "fit_pol1",
+        [](py::array_t<double, py::array::c_style | py::array::forcecast> x,
+           py::array_t<double, py::array::c_style | py::array::forcecast> y,
+           py::array_t<double, py::array::c_style | py::array::forcecast> y_err,
+           int n_threads) {
+            if (y.ndim() == 3) {
+                auto par = new NDArray<double, 3>({y.shape(0), y.shape(1), 2});
+
+                auto par_err =
+                    new NDArray<double, 3>({y.shape(0), y.shape(1), 2});
+
+                auto y_view = make_view_3d(y);
+                auto y_view_err = make_view_3d(y_err);
+                auto x_view = make_view_1d(x);
+
+                auto chi2 = new NDArray<double, 2>({y.shape(0), y.shape(1)});
+
+                aare::fit_pol1(x_view, y_view, y_view_err, par->view(),
+                               par_err->view(), chi2->view(), n_threads);
+                return py::dict("par"_a = return_image_data(par),
+                                "par_err"_a = return_image_data(par_err), 
+                                "chi2"_a = return_image_data(chi2),
+                                "Ndf"_a = y.shape(2) - 2);
+
+
+            } else if (y.ndim() == 1) {
+                auto par = new NDArray<double, 1>({2});
+                auto par_err = new NDArray<double, 1>({2});
+
+                auto y_view = make_view_1d(y);
+                auto y_view_err = make_view_1d(y_err);
+                auto x_view = make_view_1d(x);
+
+                double chi2 = 0;
+
+                aare::fit_pol1(x_view, y_view, y_view_err, par->view(),
+                               par_err->view(), chi2);
+                return py::dict("par"_a = return_image_data(par),
+                                "par_err"_a = return_image_data(par_err),
+                                "chi2"_a = chi2, "Ndf"_a = y.size() - 2);
+
+            } else {
+                throw std::runtime_error("Data must be 1D or 3D");
+            }
+        },
+        R"(
+Fit a 1D polynomial to data with error estimates.
+
+Parameters
+----------
+x : array_like
+    The x values.
+y : array_like
+    The y values.
+y_err : array_like
+    The error in the y values.
+n_threads : int, optional
+    The number of threads to use. Default is 4.
+)",
+        py::arg("x"), py::arg("y"), py::arg("y_err"), py::arg("n_threads") = 4);
+}

python/src/module.cpp

@@ -0,0 +1,35 @@
+//Files with bindings to the different classes
+#include "file.hpp"
+#include "raw_file.hpp"
+#include "ctb_raw_file.hpp"
+#include "raw_master_file.hpp"
+#include "var_cluster.hpp"
+#include "pixel_map.hpp"
+#include "pedestal.hpp"
+#include "cluster.hpp"
+#include "cluster_file.hpp"
+#include "fit.hpp"
+
+//Pybind stuff
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+namespace py = pybind11;
+
+PYBIND11_MODULE(_aare, m) {
+    define_file_io_bindings(m);
+    define_raw_file_io_bindings(m);
+    define_ctb_raw_file_io_bindings(m);
+    define_raw_master_file_bindings(m);
+    define_var_cluster_finder_bindings(m);
+    define_pixel_map_bindings(m);
+    define_pedestal_bindings<double>(m, "Pedestal_d");
+    define_pedestal_bindings<float>(m, "Pedestal_f");
+    define_cluster_finder_bindings(m);
+    define_cluster_finder_mt_bindings(m);
+    define_cluster_file_io_bindings(m);
+    define_cluster_collector_bindings(m);
+    define_cluster_file_sink_bindings(m);
+    define_fit_bindings(m);
+
+}

python/src/np_helper.hpp

@@ -0,0 +1,64 @@
+#pragma once
+
+#include <iostream>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+#include "aare/Frame.hpp"
+#include "aare/NDArray.hpp"
+#include "aare/NDView.hpp"
+
+namespace py = pybind11;
+
+// Pass image data back to python as a numpy array
+template <typename T, int64_t Ndim>
+py::array return_image_data(aare::NDArray<T, Ndim> *image) {
+
+    py::capsule free_when_done(image, [](void *f) {
+        aare::NDArray<T, Ndim> *foo =
+            reinterpret_cast<aare::NDArray<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
+}
+
+// todo rewrite generic
+template <class T, int Flags> auto get_shape_3d(py::array_t<T, Flags> arr) {
+    return aare::Shape<3>{arr.shape(0), arr.shape(1), arr.shape(2)};
+}
+
+template <class T, int Flags> auto make_view_3d(py::array_t<T, Flags> arr) {
+    return aare::NDView<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 aare::Shape<2>{arr.shape(0), arr.shape(1)};
+}
+
+template <class T, int Flags> auto get_shape_1d(py::array_t<T, Flags> arr) {
+    return aare::Shape<1>{arr.shape(0)};
+}
+
+template <class T, int Flags> auto make_view_2d(py::array_t<T, Flags> arr) {
+    return aare::NDView<T, 2>(arr.mutable_data(), get_shape_2d<T, Flags>(arr));
+}
+template <class T, int Flags> auto make_view_1d(py::array_t<T, Flags> arr) {
+    return aare::NDView<T, 1>(arr.mutable_data(), get_shape_1d<T, Flags>(arr));
+}

python/src/pedestal.hpp

@@ -0,0 +1,52 @@
+
+#include "aare/Pedestal.hpp"
+#include "np_helper.hpp"
+
+#include <cstdint>
+#include <filesystem>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+namespace py = pybind11;
+
+template <typename SUM_TYPE> void define_pedestal_bindings(py::module &m, const std::string &name) {
+    py::class_<Pedestal<SUM_TYPE>>(m, name.c_str())
+        .def(py::init<int, int, int>())
+        .def(py::init<int, int>())
+        .def("mean",
+             [](Pedestal<SUM_TYPE> &self) {
+                 auto mea = new NDArray<SUM_TYPE, 2>{};
+                 *mea = self.mean();
+                 return return_image_data(mea);
+             })
+        .def("variance", [](Pedestal<SUM_TYPE> &self) {
+            auto var = new NDArray<SUM_TYPE, 2>{};
+            *var = self.variance();
+            return return_image_data(var);
+        })
+        .def("std", [](Pedestal<SUM_TYPE> &self) {
+            auto std = new NDArray<SUM_TYPE, 2>{};
+            *std = self.std();
+            return return_image_data(std);
+        })
+        .def("clear", py::overload_cast<>(&Pedestal<SUM_TYPE>::clear))
+        .def_property_readonly("rows", &Pedestal<SUM_TYPE>::rows)
+        .def_property_readonly("cols", &Pedestal<SUM_TYPE>::cols)
+        .def_property_readonly("n_samples", &Pedestal<SUM_TYPE>::n_samples)
+        .def_property_readonly("sum", &Pedestal<SUM_TYPE>::get_sum)
+        .def_property_readonly("sum2", &Pedestal<SUM_TYPE>::get_sum2)
+        .def("clone",
+             [&](Pedestal<SUM_TYPE> &pedestal) {
+                 return Pedestal<SUM_TYPE>(pedestal);
+             })
+        //TODO! add push for other data types
+        .def("push", [](Pedestal<SUM_TYPE> &pedestal, py::array_t<uint16_t> &f) {
+            auto v = make_view_2d(f);
+            pedestal.push(v);
+        })
+        .def("push_no_update", [](Pedestal<SUM_TYPE> &pedestal, py::array_t<uint16_t, py::array::c_style> &f) {
+            auto v = make_view_2d(f);
+            pedestal.push_no_update(v);
+        }, py::arg().noconvert())
+        .def("update_mean", &Pedestal<SUM_TYPE>::update_mean);
+}

python/src/pixel_map.hpp

@@ -0,0 +1,41 @@
+#include "aare/PixelMap.hpp"
+#include "np_helper.hpp"
+
+
+#include <cstdint>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+
+
+namespace py = pybind11;
+using namespace::aare;
+
+
+void define_pixel_map_bindings(py::module &m) {
+    m.def("GenerateMoench03PixelMap", []() {
+        auto ptr = new NDArray<ssize_t,2>(GenerateMoench03PixelMap());
+        return return_image_data(ptr);
+    })
+    .def("GenerateMoench05PixelMap", []() {
+        auto ptr = new NDArray<ssize_t,2>(GenerateMoench05PixelMap());
+        return return_image_data(ptr);
+    })
+        .def("GenerateMoench05PixelMap1g", []() {
+        auto ptr = new NDArray<ssize_t,2>(GenerateMoench05PixelMap1g());
+        return return_image_data(ptr);
+    })
+    .def("GenerateMoench05PixelMapOld", []() {
+        auto ptr = new NDArray<ssize_t,2>(GenerateMoench05PixelMapOld());
+        return return_image_data(ptr);
+    })
+    .def("GenerateMH02SingleCounterPixelMap", []() {
+        auto ptr = new NDArray<ssize_t,2>(GenerateMH02SingleCounterPixelMap());
+        return return_image_data(ptr);
+    })
+    .def("GenerateMH02FourCounterPixelMap", []() {
+        auto ptr = new NDArray<ssize_t,3>(GenerateMH02FourCounterPixelMap());
+        return return_image_data(ptr);
+    });
+
+}

python/src/raw_file.hpp

@@ -0,0 +1,106 @@
+#include "aare/CtbRawFile.hpp"
+#include "aare/File.hpp"
+#include "aare/Frame.hpp"
+#include "aare/RawFile.hpp"
+#include "aare/RawMasterFile.hpp"
+#include "aare/RawSubFile.hpp"
+
+#include "aare/defs.hpp"
+// #include "aare/fClusterFileV2.hpp"
+
+#include <cstdint>
+#include <filesystem>
+#include <pybind11/iostream.h>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl/filesystem.h>
+#include <string>
+
+namespace py = pybind11;
+using namespace ::aare;
+
+void define_raw_file_io_bindings(py::module &m) {
+    py::class_<RawFile>(m, "RawFile")
+        .def(py::init<const std::filesystem::path &>())
+        .def("read_frame",
+             [](RawFile &self) {
+                 py::array image;
+                 std::vector<ssize_t> shape;
+                 shape.reserve(2);
+                 shape.push_back(self.rows());
+                 shape.push_back(self.cols());
+
+                 // return headers from all subfiles
+                 py::array_t<DetectorHeader> header(self.n_mod());
+
+                 const uint8_t item_size = self.bytes_per_pixel();
+                 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()),
+                     header.mutable_data());
+
+                 return py::make_tuple(header, image);
+             })
+        .def(
+            "read_n",
+            [](RawFile &self, size_t n_frames) {
+                // adjust for actual frames left in the file
+                n_frames =
+                    std::min(n_frames, self.total_frames() - self.tell());
+                if (n_frames == 0) {
+                    throw std::runtime_error("No frames left in file");
+                }
+                std::vector<size_t> shape{n_frames, self.rows(), self.cols()};
+            
+                // return headers from all subfiles
+                py::array_t<DetectorHeader> header;
+                if (self.n_mod() == 1) {
+                    header = py::array_t<DetectorHeader>(n_frames);
+                } else {
+                    header = py::array_t<DetectorHeader>({self.n_mod(), n_frames});
+                }
+                // py::array_t<DetectorHeader> header({self.n_mod(), n_frames});
+
+                py::array image;
+                const uint8_t item_size = self.bytes_per_pixel();
+                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()),
+                    n_frames, header.mutable_data());
+
+                return py::make_tuple(header, image);
+            },
+            R"(
+             Read n frames from the file.
+             )")
+        .def("frame_number", &RawFile::frame_number)
+        .def_property_readonly("bytes_per_frame", &RawFile::bytes_per_frame)
+        .def_property_readonly("pixels_per_frame", &RawFile::pixels_per_frame)
+        .def_property_readonly("bytes_per_pixel", &RawFile::bytes_per_pixel)
+        .def("seek", &RawFile::seek, R"(
+            Seek to a frame index in file.
+            )")
+        .def("tell", &RawFile::tell, R"(
+            Return the current frame number.)")
+        .def_property_readonly("total_frames", &RawFile::total_frames)
+        .def_property_readonly("rows", &RawFile::rows)
+        .def_property_readonly("cols", &RawFile::cols)
+        .def_property_readonly("bitdepth", &RawFile::bitdepth)
+        .def_property_readonly("geometry", &RawFile::geometry)
+        .def_property_readonly("n_mod", &RawFile::n_mod)
+        .def_property_readonly("detector_type", &RawFile::detector_type)
+        .def_property_readonly("master", &RawFile::master);
+}

python/src/raw_master_file.hpp

@@ -0,0 +1,85 @@
+
+#include "aare/CtbRawFile.hpp"
+#include "aare/File.hpp"
+#include "aare/Frame.hpp"
+#include "aare/RawFile.hpp"
+#include "aare/RawMasterFile.hpp"
+#include "aare/RawSubFile.hpp"
+
+#include "aare/defs.hpp"
+// #include "aare/fClusterFileV2.hpp"
+
+#include <cstdint>
+#include <filesystem>
+#include <pybind11/iostream.h>
+#include <pybind11/numpy.h>
+#include <pybind11/pybind11.h>
+#include <pybind11/stl.h>
+#include <pybind11/stl/filesystem.h>
+#include <string>
+
+namespace py = pybind11;
+using namespace ::aare;
+
+void define_raw_master_file_bindings(py::module &m) {
+    py::class_<RawMasterFile>(m, "RawMasterFile")
+        .def(py::init<const std::filesystem::path &>())
+        .def("data_fname", &RawMasterFile::data_fname, R"(
+
+            Parameters
+            ------------
+                module_index : int 
+                        module index (d0, d1 .. dN)
+                file_index : int 
+                        file index (f0, f1 .. fN)
+
+            Returns
+            ----------
+                os.PathLike
+                        The name of the data file.
+
+            )")
+        .def_property_readonly("version", &RawMasterFile::version)
+        .def_property_readonly("detector_type", &RawMasterFile::detector_type)
+        .def_property_readonly("timing_mode", &RawMasterFile::timing_mode)
+        .def_property_readonly("image_size_in_bytes",
+                               &RawMasterFile::image_size_in_bytes)
+        .def_property_readonly("frames_in_file", &RawMasterFile::frames_in_file)
+        .def_property_readonly("pixels_y", &RawMasterFile::pixels_y)
+        .def_property_readonly("pixels_x", &RawMasterFile::pixels_x)
+        .def_property_readonly("max_frames_per_file",
+                               &RawMasterFile::max_frames_per_file)
+        .def_property_readonly("bitdepth", &RawMasterFile::bitdepth)
+        .def_property_readonly("frame_padding", &RawMasterFile::frame_padding)
+        .def_property_readonly("frame_discard_policy",
+                               &RawMasterFile::frame_discard_policy)
+
+        .def_property_readonly("total_frames_expected",
+                               &RawMasterFile::total_frames_expected)
+        .def_property_readonly("geometry", &RawMasterFile::geometry)
+        .def_property_readonly("analog_samples", &RawMasterFile::analog_samples, R"(
+            Number of analog samples
+
+            Returns
+            ----------
+                int | None
+                    The number of analog samples in the file (or None if not enabled)
+            )")
+        .def_property_readonly("digital_samples",
+                               &RawMasterFile::digital_samples,  R"(
+            Number of digital samples
+
+            Returns
+            ----------
+                int | None
+                    The number of digital samples in the file (or None if not enabled)
+            )")
+
+        .def_property_readonly("transceiver_samples",
+                               &RawMasterFile::transceiver_samples)
+        .def_property_readonly("number_of_rows", &RawMasterFile::number_of_rows)
+        .def_property_readonly("quad", &RawMasterFile::quad)
+        .def_property_readonly("scan_parameters",
+                               &RawMasterFile::scan_parameters)
+        .def_property_readonly("roi", &RawMasterFile::roi);
+}