fixed calculate eta

Fixed warnings:
- unused variable in Interpolator
- Narrowing conversions uint64-->int64

Removed an ambiguous function from JungfrauDataFile
- NDarry read_frame(header&=nullptr)
- Frame read_frame()

NDArray and NDView size() is now signed

.clang-format

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

.clang-tidy

@@ -0,0 +1,42 @@
+
+---
+Checks:          '*, 
+                    -altera-*,
+                    -android-cloexec-fopen,
+                    -cppcoreguidelines-pro-bounds-array-to-pointer-decay, 
+                    -cppcoreguidelines-pro-bounds-pointer-arithmetic, 
+                    -fuchsia*,
+                    -readability-else-after-return,
+                    -readability-avoid-const-params-in-decls,
+                    -readability-identifier-length,
+                    -cppcoreguidelines-pro-bounds-constant-array-index,
+                    -cppcoreguidelines-pro-type-reinterpret-cast,
+                    -llvm-header-guard,
+                    -modernize-use-nodiscard,
+                    -misc-non-private-member-variables-in-classes,
+                    -readability-static-accessed-through-instance,
+                    -readability-braces-around-statements,
+                    -readability-isolate-declaration,
+                    -readability-implicit-bool-conversion,
+                    -readability-identifier-length,
+                    -readability-identifier-naming,
+                    -hicpp-signed-bitwise,
+                    -hicpp-no-array-decay, 
+                    -hicpp-braces-around-statements,
+                    -google-runtime-references, 
+                    -google-readability-todo, 
+                    -google-readability-braces-around-statements,
+                    -modernize-use-trailing-return-type,
+                    -llvmlibc-*'
+
+HeaderFilterRegex: \.hpp
+FormatStyle:     none
+CheckOptions:
+  - { key: readability-identifier-naming.NamespaceCase,       value: lower_case }
+ # - { key: readability-identifier-naming.FunctionCase,        value: lower_case }
+  - { key: readability-identifier-naming.ClassCase,           value: CamelCase  }
+ # - { key: readability-identifier-naming.MethodCase,          value: CamelCase  }
+ # - { key: readability-identifier-naming.StructCase,          value: CamelCase  }
+ # - { key: readability-identifier-naming.VariableCase,        value: lower_case }
+  - { key: readability-identifier-naming.GlobalConstantCase,  value: UPPER_CASE }
+...

.gitea/workflows/cmake_build.yml

@@ -0,0 +1,58 @@
+name: Build the package using cmake then documentation
+
+on:
+  workflow_dispatch:
+
+    
+
+permissions:
+  contents: read
+  pages: write
+  id-token: write
+
+jobs:
+  build:
+    strategy:
+      fail-fast: false
+      matrix:
+        platform: [ubuntu-latest, ] 
+        python-version: ["3.12", ] 
+
+    runs-on: ${{ matrix.platform }}
+
+
+    defaults:
+      run:
+        shell: "bash -l {0}"
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Setup dev env
+        run: |
+            sudo apt-get update
+            sudo apt-get -y install cmake gcc g++
+
+      - name: Get conda
+        uses: conda-incubator/setup-miniconda@v3
+        with:
+          python-version: ${{ matrix.python-version }}
+          environment-file: etc/dev-env.yml
+          miniforge-version: latest
+          channels: conda-forge
+          conda-remove-defaults: "true"
+
+      - name: Build library 
+        run: |
+            mkdir build 
+            cd build 
+            cmake .. -DAARE_SYSTEM_LIBRARIES=ON -DAARE_DOCS=ON 
+            make -j 2 
+            make docs
+
+
+
+
+
+
+

.gitea/workflows/rh8-native.yml

@@ -0,0 +1,30 @@
+name: Build on RHEL8
+
+on:
+  workflow_dispatch:
+
+permissions:
+  contents: read
+
+jobs:
+  buildh:
+    runs-on: "ubuntu-latest"
+    container:
+      image: gitea.psi.ch/images/rhel8-developer-gitea-actions
+    steps:
+      - uses: actions/checkout@v4
+
+
+      - name: Install dependencies
+        run: |
+          dnf install -y cmake python3.12 python3.12-devel python3.12-pip
+
+      - name: Build library 
+        run: |
+            mkdir build && cd build
+            cmake .. -DAARE_PYTHON_BINDINGS=ON -DAARE_TESTS=ON
+            make -j 2 
+
+      - name: C++ unit tests
+        working-directory: ${{gitea.workspace}}/build
+        run: ctest

.gitea/workflows/rh9-native.yml

@@ -0,0 +1,31 @@
+name: Build on RHEL9
+
+on:
+  push:
+  workflow_dispatch:
+
+permissions:
+  contents: read
+
+jobs:
+  buildh:
+    runs-on: "ubuntu-latest"
+    container:
+      image: gitea.psi.ch/images/rhel9-developer-gitea-actions
+    steps:
+      - uses: actions/checkout@v4
+
+
+      - name: Install dependencies
+        run: |
+          dnf install -y cmake python3.12 python3.12-devel python3.12-pip
+
+      - name: Build library 
+        run: |
+            mkdir build && cd build
+            cmake .. -DAARE_PYTHON_BINDINGS=ON -DAARE_TESTS=ON
+            make -j 2 
+
+      - name: C++ unit tests
+        working-directory: ${{gitea.workspace}}/build
+        run: ctest

.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,66 @@
+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, ] 
+        python-version: ["3.12",]
+
+    runs-on: ${{ matrix.platform }}
+
+    defaults:
+      run:
+        shell: "bash -l {0}"
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Get conda
+        uses: conda-incubator/setup-miniconda@v3
+        with:
+          python-version: ${{ matrix.python-version }}
+          environment-file: etc/dev-env.yml
+          miniforge-version: latest
+          channels: conda-forge
+          conda-remove-defaults: "true"
+
+      - 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

@@ -21,28 +21,37 @@ 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})
 
-option(AARE_USE_WARNINGS "Enable warnings" ON)
-option(AARE_PYTHON_BINDINGS "Build python bindings" ON)
-option(AARE_TESTS "Build tests" ON)
-option(AARE_EXAMPLES "Build examples" ON)
+
+# General options
+option(AARE_PYTHON_BINDINGS "Build python bindings" OFF)
+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(ENABLE_DRAFTS "Enable zmq drafts (depends on gnutls or nss)" OFF)
+option(AARE_FETCH_LMFIT "Use FetchContent to download lmfit" ON)
 
 
-#Convenience option to use system libraries
+#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")
@@ -51,17 +60,86 @@ if(AARE_SYSTEM_LIBRARIES)
     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)
+    # Still fetch lmfit when setting AARE_SYSTEM_LIBRARIES since this is not available
+    # on conda-forge
+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)
+    #TODO! Should we fetch lmfit from the web or inlcude a tar.gz in the repo?
+    set(LMFIT_PATCH_COMMAND git apply ${CMAKE_CURRENT_SOURCE_DIR}/patches/lmfit.patch)
+
+    # For cmake < 3.28 we can't supply EXCLUDE_FROM_ALL to FetchContent_Declare
+    # so we need this workaround
+    if (${CMAKE_VERSION} VERSION_LESS "3.28")
+        FetchContent_Declare(
+            lmfit
+            GIT_REPOSITORY https://jugit.fz-juelich.de/mlz/lmfit.git
+            GIT_TAG        main
+            PATCH_COMMAND ${LMFIT_PATCH_COMMAND}
+            UPDATE_DISCONNECTED 1
+        )
+    else()
+        FetchContent_Declare(
+            lmfit
+            GIT_REPOSITORY https://jugit.fz-juelich.de/mlz/lmfit.git
+            GIT_TAG        main
+            PATCH_COMMAND ${LMFIT_PATCH_COMMAND}
+            UPDATE_DISCONNECTED 1
+            EXCLUDE_FROM_ALL 1
+        )
+    endif()
+    
+
+    #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 "")
+
+    if (${CMAKE_VERSION} VERSION_LESS "3.28")
+        if(NOT lmfit_POPULATED)
+            FetchContent_Populate(lmfit)
+            add_subdirectory(${lmfit_SOURCE_DIR} ${lmfit_BINARY_DIR} EXCLUDE_FROM_ALL)
+        endif()
+    else()
+        FetchContent_MakeAvailable(lmfit)
+    endif()
+    
+    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()
+    set(ZMQ_PATCH_COMMAND git apply ${CMAKE_CURRENT_SOURCE_DIR}/patches/libzmq_cmake_version.patch)
     FetchContent_Declare(
         libzmq
         GIT_REPOSITORY https://github.com/zeromq/libzmq.git
         GIT_TAG        v4.3.4
+        PATCH_COMMAND ${ZMQ_PATCH_COMMAND}
+        UPDATE_DISCONNECTED 1
     )
     # Disable unwanted options from libzmq
     set(BUILD_TESTS OFF CACHE BOOL "Switch off libzmq test build")
@@ -94,15 +172,39 @@ if (AARE_FETCH_FMT)
             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()
 
-#TODO! Add options for nlohmann_json as well
-find_package(nlohmann_json 3.11.3 REQUIRED)
 
+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)
 
@@ -176,42 +278,45 @@ if(CMAKE_BUILD_TYPE STREQUAL "Release")
     target_compile_options(aare_compiler_flags INTERFACE -O3)
 else()
     message(STATUS "Debug build")
-    target_compile_options(
-        aare_compiler_flags 
-        INTERFACE 
-            -Og 
-            -ggdb3 
-            # -D_GLIBCXX_DEBUG # causes errors with boost
-            -D_GLIBCXX_DEBUG_PEDANTIC
+endif()
+
+# 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()
-
-if(AARE_USE_WARNINGS)
-        target_compile_options(
-            aare_compiler_flags 
-            INTERFACE 
-                -Wall 
-                -Wextra 
-                -pedantic 
-                -Wshadow 
-                -Wnon-virtual-dtor
-                -Woverloaded-virtual
-                -Wdouble-promotion
-                -Wformat=2
-                -Wredundant-decls
-                -Wvla
-                -Wdouble-promotion
-                -Werror=return-type #important can cause segfault in optimzed builds
-            )
-endif()
-
-
-
 endif() #GCC/Clang specific
 
 
-
+if(AARE_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()
 
 
 
@@ -226,47 +331,79 @@ endif()
 ###------------------------------------------------------------------------------------------
 
 set(PUBLICHEADERS
+    include/aare/ArrayExpr.hpp
+    include/aare/CalculateEta.hpp
+    include/aare/Cluster.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/FilePtr.hpp
     include/aare/Frame.hpp
-    include/aare/json.hpp
+    include/aare/GainMap.hpp
+    include/aare/geo_helpers.hpp
+    include/aare/JungfrauDataFile.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/SubFile.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/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/FilePtr.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/Fit.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/geo_helpers.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/JungfrauDataFile.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.cpp
-    ${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.cpp
-    ${CMAKE_CURRENT_SOURCE_DIR}/src/SubFile.cpp
     ${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/Interpolator.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/PixelMap.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/RawSubFile.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.cpp
+    ${CMAKE_CURRENT_SOURCE_DIR}/src/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}>"
+    "$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>" 
 )
 
-target_link_libraries(aare_core PUBLIC fmt::fmt PRIVATE aare_compiler_flags nlohmann_json::nlohmann_json)
+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_NAME SlsDetectorStatic
     ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}
     PUBLIC_HEADER "${PUBLICHEADERS}"
 )
@@ -277,16 +414,25 @@ endif()
 
 if(AARE_TESTS)
     set(TestSources
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/algorithm.test.cpp
             ${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/Cluster.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/CalculateEta.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFile.test.cpp
             ${CMAKE_CURRENT_SOURCE_DIR}/src/Pedestal.test.cpp
+            ${CMAKE_CURRENT_SOURCE_DIR}/src/JungfrauDataFile.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} )
@@ -312,13 +458,13 @@ 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>
-)
+# #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)
 
@@ -354,7 +500,7 @@ 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 {}"
+    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
@@ -362,6 +508,6 @@ add_custom_target(
 
 if(AARE_MASTER_PROJECT)
     set(CMAKE_INSTALL_DIR "share/cmake/${PROJECT_NAME}")
-    set(PROJECT_LIBRARIES slsSupportShared slsDetectorShared slsReceiverShared)
+    set(PROJECT_LIBRARIES aare-core aare-compiler-flags )
     include(cmake/package_config.cmake)
-endif()
+endif()

README.md

@@ -2,19 +2,70 @@
 Data analysis library for PSI hybrid detectors
 
 
-## Status
+## Build and install
 
-- [ ] Build with CMake on RH8
-- [ ] conda package
+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
 
 
-## Project structure 
+### Install using conda/mamba
 
-include/aare - public headers
+```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
 
 
-## Open questions
+#Now configure your project
+ cmake .. -DCMAKE_PREFIX_PATH=SOME_PATH
+```
 
-- How many sub libraries? 
-- Where to place test data? This data is also needed for github actions...
-- What to return to numpy? Our NDArray or a numpy ndarray? Lifetime? 
+### Local build of conda pkgs
+
+```bash
+conda build . --variants="{python: [3.11, 3.12, 3.13]}"
+```

benchmarks/CMakeLists.txt

@@ -0,0 +1,27 @@
+
+include(FetchContent)
+
+
+FetchContent_Declare(
+    benchmark
+    GIT_REPOSITORY https://github.com/google/benchmark.git
+    GIT_TAG        v1.8.3 # Change to the latest version if needed
+)
+
+# Ensure Google Benchmark is built correctly
+set(BENCHMARK_ENABLE_TESTING OFF CACHE BOOL "" FORCE)
+
+FetchContent_MakeAvailable(benchmark)
+
+add_executable(benchmarks)
+
+target_sources(benchmarks PRIVATE ndarray_benchmark.cpp calculateeta_benchmark.cpp)
+
+# Link Google Benchmark and other necessary libraries
+target_link_libraries(benchmarks PRIVATE benchmark::benchmark aare_core aare_compiler_flags)
+
+# Set output properties
+set_target_properties(benchmarks PROPERTIES
+    RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}
+    OUTPUT_NAME run_benchmarks
+)

benchmarks/calculateeta_benchmark.cpp

@@ -0,0 +1,70 @@
+#include "aare/CalculateEta.hpp"
+#include "aare/ClusterFile.hpp"
+#include <benchmark/benchmark.h>
+
+using namespace aare;
+
+class ClusterFixture : public benchmark::Fixture {
+  public:
+    Cluster<int, 2, 2> cluster_2x2{};
+    Cluster<int, 3, 3> cluster_3x3{};
+
+  private:
+    using benchmark::Fixture::SetUp;
+
+    void SetUp([[maybe_unused]] const benchmark::State &state) override {
+        int temp_data[4] = {1, 2, 3, 1};
+        std::copy(std::begin(temp_data), std::end(temp_data),
+                  std::begin(cluster_2x2.data));
+
+        cluster_2x2.x = 0;
+        cluster_2x2.y = 0;
+
+        int temp_data2[9] = {1, 2, 3, 1, 3, 4, 5, 1, 20};
+        std::copy(std::begin(temp_data2), std::end(temp_data2),
+                  std::begin(cluster_3x3.data));
+
+        cluster_3x3.x = 0;
+        cluster_3x3.y = 0;
+    }
+
+    // void TearDown(::benchmark::State& state) {
+    // }
+};
+
+BENCHMARK_F(ClusterFixture, Calculate2x2Eta)(benchmark::State &st) {
+    for (auto _ : st) {
+        // This code gets timed
+        Eta2 eta = calculate_eta2(cluster_2x2);
+        benchmark::DoNotOptimize(eta);
+    }
+}
+
+// almost takes double the time
+BENCHMARK_F(ClusterFixture,
+            CalculateGeneralEtaFor2x2Cluster)(benchmark::State &st) {
+    for (auto _ : st) {
+        // This code gets timed
+        Eta2 eta = calculate_eta2<int, 2, 2>(cluster_2x2);
+        benchmark::DoNotOptimize(eta);
+    }
+}
+
+BENCHMARK_F(ClusterFixture, Calculate3x3Eta)(benchmark::State &st) {
+    for (auto _ : st) {
+        // This code gets timed
+        Eta2 eta = calculate_eta2(cluster_3x3);
+        benchmark::DoNotOptimize(eta);
+    }
+}
+
+// almost takes double the time
+BENCHMARK_F(ClusterFixture,
+            CalculateGeneralEtaFor3x3Cluster)(benchmark::State &st) {
+    for (auto _ : st) {
+        // This code gets timed
+        Eta2 eta = calculate_eta2<int, 3, 3>(cluster_3x3);
+        benchmark::DoNotOptimize(eta);
+    }
+}
+// BENCHMARK_MAIN();

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/helpers.cmake

@@ -4,3 +4,43 @@ if (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
   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/project-config.cmake.in

@@ -12,8 +12,10 @@ 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)

conda-recipe/build.sh

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

conda-recipe/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/copy_lib.sh

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

conda-recipe/meta.yaml

@@ -1,104 +1,55 @@
-
-
 package:
-  name: aare_software
-  version: {{ environ.get('GIT_DESCRIBE_TAG', '') }}
+  name: aare
+  version: 2025.4.1 #TODO! how to not duplicate this?
+
+
+
+
 
 source:
-  - path: ..
+  path: ..
 
 build:
   number: 0
-  binary_relocation: True
-  rpaths:
-    - lib/
+  script:
+    - unset CMAKE_GENERATOR && {{ PYTHON }} -m pip install . -vv  # [not win]
+    - {{ PYTHON }} -m pip install . -vv  # [win]
 
 requirements:
   build:
-    - {{ compiler('c') }}
-    - {{compiler('cxx')}}
-    - cmake
-    # - qt 5.*
-    # - xorg-libx11
-    # - xorg-libice
-    # - xorg-libxext
-    # - xorg-libsm
-    # - xorg-libxau
-    # - xorg-libxrender
-    # - xorg-libxfixes
-    # - {{ cdt('mesa-libgl-devel') }}      # [linux]
-    # - {{ cdt('mesa-libegl-devel') }}     # [linux]
-    # - {{ cdt('mesa-dri-drivers') }}      # [linux]
-    # - {{ cdt('libselinux') }}            # [linux]
-    # - {{ cdt('libxdamage') }}            # [linux]
-    # - {{ cdt('libxxf86vm') }}            # [linux]
-    # - expat
+    - python {{python}}
+    - numpy {{ numpy }}
+    - {{ compiler('cxx') }}
+
 
   host:
-    # - libstdcxx-ng
-    # - libgcc-ng
-    # - xorg-libx11
-    # - xorg-libice
-    # - xorg-libxext
-    # - xorg-libsm
-    # - xorg-libxau
-    # - xorg-libxrender
-    # - xorg-libxfixes
-    # - expat
+    - cmake
+    - ninja
+    - python {{python}}
+    - numpy {{ numpy }}
+    - pip
+    - scikit-build-core
+    - pybind11 >=2.13.0
+    - fmt
+    - zeromq
+    - nlohmann_json
+    - catch2
 
   run:
-    # - libstdcxx-ng
-    # - libgcc-ng 
+    - python {{python}}
+    - numpy {{ numpy }}
 
 
-outputs:
-  - name: aarelib
-    script: copy_lib.sh
+test:
+  imports:
+    - aare
+  # requires:
+  #   - pytest
+  # source_files:
+  #   - tests
+  # commands:
+  #   - pytest tests
 
-    requirements:
-      build:
-        - {{ compiler('c') }}
-        - {{compiler('cxx')}}
-        - catch2
-        - zstd
-
-        # - libstdcxx-ng
-        # - libgcc-ng
-
-
-      run: 
-        # - libstdcxx-ng
-        # - libgcc-ng
-
-  # - name: aare
-
-  #   script: build_pylib.sh
-
-  #   requirements:
-  #     build:
-  #       - python
-  #       - {{ compiler('c') }}
-  #       - {{compiler('cxx')}}
-  #       - {{ pin_subpackage('slsdetlib', exact=True) }}
-  #       - setuptools
-  #       - pybind11=2.11
-
-  #     host:
-  #       - python
-  #       - {{ pin_subpackage('slsdetlib', exact=True) }}
-  #       - pybind11=2.11
-
-
-  #     run:
-  #       - libstdcxx-ng
-  #       - libgcc-ng
-  #       - python
-  #       - numpy
-  #       - {{ pin_subpackage('slsdetlib', exact=True) }}
-
-
-  #   test:
-  #     imports:
-  #       - slsdet
-
-  
+about:
+  summary: An example project built with pybind11 and scikit-build.
+  # license_file: LICENSE

docs/CMakeLists.txt

@@ -10,22 +10,15 @@ configure_file(${DOXYGEN_IN} ${DOXYGEN_OUT} @ONLY)
 set(SPHINX_SOURCE ${CMAKE_CURRENT_SOURCE_DIR}/src)
 set(SPHINX_BUILD ${CMAKE_CURRENT_BINARY_DIR})
 
-set(SPHINX_SOURCE_FILES
-  src/index.rst
-  src/NDArray.rst
-  src/NDView.rst
-  src/File.rst
-  src/Frame.rst
-  src/Dtype.rst
-  src/ClusterFinder.rst
-  src/Pedestal.rst
-  src/VarClusterFinder.rst
-  src/pyFile.rst
-)
+
+file(GLOB SPHINX_SOURCE_FILES CONFIGURE_DEPENDS "src/*.rst")
+
+
 
 foreach(filename ${SPHINX_SOURCE_FILES})
-configure_file(${CMAKE_CURRENT_SOURCE_DIR}/${filename} 
-              "${SPHINX_BUILD}/${filename}")
+    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(

docs/Doxyfile.in

@@ -241,11 +241,7 @@ TAB_SIZE               = 4
 
 ALIASES                =
 
-# This tag can be used to specify a number of word-keyword mappings (TCL only).
-# A mapping has the form "name=value". For example adding "class=itcl::class"
-# will allow you to use the command class in the itcl::class meaning.
 
-TCL_SUBST              =
 
 # Set the OPTIMIZE_OUTPUT_FOR_C tag to YES if your project consists of C sources
 # only. Doxygen will then generate output that is more tailored for C. For
@@ -1082,12 +1078,7 @@ VERBATIM_HEADERS       = YES
 
 ALPHABETICAL_INDEX     = YES
 
-# The COLS_IN_ALPHA_INDEX tag can be used to specify the number of columns in
-# which the alphabetical index list will be split.
-# Minimum value: 1, maximum value: 20, default value: 5.
-# This tag requires that the tag ALPHABETICAL_INDEX is set to YES.
 
-COLS_IN_ALPHA_INDEX    = 5
 
 # In case all classes in a project start with a common prefix, all classes will
 # be put under the same header in the alphabetical index. The IGNORE_PREFIX tag
@@ -1216,14 +1207,6 @@ HTML_COLORSTYLE_SAT    = 100
 
 HTML_COLORSTYLE_GAMMA  = 80
 
-# If the HTML_TIMESTAMP tag is set to YES then the footer of each generated HTML
-# page will contain the date and time when the page was generated. Setting this
-# to YES can help to show when doxygen was last run and thus if the
-# documentation is up to date.
-# The default value is: NO.
-# This tag requires that the tag GENERATE_HTML is set to YES.
-
-HTML_TIMESTAMP         = NO
 
 # If the HTML_DYNAMIC_MENUS tag is set to YES then the generated HTML
 # documentation will contain a main index with vertical navigation menus that
@@ -1248,7 +1231,7 @@ HTML_DYNAMIC_SECTIONS  = NO
 # shown in the various tree structured indices initially; the user can expand
 # and collapse entries dynamically later on. Doxygen will expand the tree to
 # such a level that at most the specified number of entries are visible (unless
-# a fully collapsed tree already exceeds this amount). So setting the number of
+# a fully collapsed tree exceeds this amount). So setting the number of
 # entries 1 will produce a full collapsed tree by default. 0 is a special value
 # representing an infinite number of entries and will result in a full expanded
 # tree by default.
@@ -1503,16 +1486,6 @@ EXT_LINKS_IN_WINDOW    = NO
 
 FORMULA_FONTSIZE       = 10
 
-# Use the FORMULA_TRANSPARENT tag to determine whether or not the images
-# generated for formulas are transparent PNGs. Transparent PNGs are not
-# supported properly for IE 6.0, but are supported on all modern browsers.
-#
-# Note that when changing this option you need to delete any form_*.png files in
-# the HTML output directory before the changes have effect.
-# The default value is: YES.
-# This tag requires that the tag GENERATE_HTML is set to YES.
-
-FORMULA_TRANSPARENT    = YES
 
 # Enable the USE_MATHJAX option to render LaTeX formulas using MathJax (see
 # https://www.mathjax.org) which uses client side Javascript for the rendering
@@ -1776,31 +1749,6 @@ PDF_HYPERLINKS         = YES
 
 USE_PDFLATEX           = YES
 
-# If the LATEX_BATCHMODE tag is set to YES, doxygen will add the \batchmode
-# command to the generated LaTeX files. This will instruct LaTeX to keep running
-# if errors occur, instead of asking the user for help. This option is also used
-# when generating formulas in HTML.
-# The default value is: NO.
-# This tag requires that the tag GENERATE_LATEX is set to YES.
-
-LATEX_BATCHMODE        = NO
-
-# If the LATEX_HIDE_INDICES tag is set to YES then doxygen will not include the
-# index chapters (such as File Index, Compound Index, etc.) in the output.
-# The default value is: NO.
-# This tag requires that the tag GENERATE_LATEX is set to YES.
-
-LATEX_HIDE_INDICES     = NO
-
-# If the LATEX_SOURCE_CODE tag is set to YES then doxygen will include source
-# code with syntax highlighting in the LaTeX output.
-#
-# Note that which sources are shown also depends on other settings such as
-# SOURCE_BROWSER.
-# The default value is: NO.
-# This tag requires that the tag GENERATE_LATEX is set to YES.
-
-LATEX_SOURCE_CODE      = NO
 
 # The LATEX_BIB_STYLE tag can be used to specify the style to use for the
 # bibliography, e.g. plainnat, or ieeetr. See
@@ -1810,227 +1758,6 @@ LATEX_SOURCE_CODE      = NO
 
 LATEX_BIB_STYLE        = plain
 
-# If the LATEX_TIMESTAMP tag is set to YES then the footer of each generated
-# page will contain the date and time when the page was generated. Setting this
-# to NO can help when comparing the output of multiple runs.
-# The default value is: NO.
-# This tag requires that the tag GENERATE_LATEX is set to YES.
-
-LATEX_TIMESTAMP        = NO
-
-#---------------------------------------------------------------------------
-# Configuration options related to the RTF output
-#---------------------------------------------------------------------------
-
-# If the GENERATE_RTF tag is set to YES, doxygen will generate RTF output. The
-# RTF output is optimized for Word 97 and may not look too pretty with other RTF
-# readers/editors.
-# The default value is: NO.
-
-GENERATE_RTF           = NO
-
-# The RTF_OUTPUT tag is used to specify where the RTF docs will be put. If a
-# relative path is entered the value of OUTPUT_DIRECTORY will be put in front of
-# it.
-# The default directory is: rtf.
-# This tag requires that the tag GENERATE_RTF is set to YES.
-
-RTF_OUTPUT             = rtf
-
-# If the COMPACT_RTF tag is set to YES, doxygen generates more compact RTF
-# documents. This may be useful for small projects and may help to save some
-# trees in general.
-# The default value is: NO.
-# This tag requires that the tag GENERATE_RTF is set to YES.
-
-COMPACT_RTF            = NO
-
-# If the RTF_HYPERLINKS tag is set to YES, the RTF that is generated will
-# contain hyperlink fields. The RTF file will contain links (just like the HTML
-# output) instead of page references. This makes the output suitable for online
-# browsing using Word or some other Word compatible readers that support those
-# fields.
-#
-# Note: WordPad (write) and others do not support links.
-# The default value is: NO.
-# This tag requires that the tag GENERATE_RTF is set to YES.
-
-RTF_HYPERLINKS         = NO
-
-# Load stylesheet definitions from file. Syntax is similar to doxygen's config
-# file, i.e. a series of assignments. You only have to provide replacements,
-# missing definitions are set to their default value.
-#
-# See also section "Doxygen usage" for information on how to generate the
-# default style sheet that doxygen normally uses.
-# This tag requires that the tag GENERATE_RTF is set to YES.
-
-RTF_STYLESHEET_FILE    =
-
-# Set optional variables used in the generation of an RTF document. Syntax is
-# similar to doxygen's config file. A template extensions file can be generated
-# using doxygen -e rtf extensionFile.
-# This tag requires that the tag GENERATE_RTF is set to YES.
-
-RTF_EXTENSIONS_FILE    =
-
-# If the RTF_SOURCE_CODE tag is set to YES then doxygen will include source code
-# with syntax highlighting in the RTF output.
-#
-# Note that which sources are shown also depends on other settings such as
-# SOURCE_BROWSER.
-# The default value is: NO.
-# This tag requires that the tag GENERATE_RTF is set to YES.
-
-RTF_SOURCE_CODE        = NO
-
-#---------------------------------------------------------------------------
-# Configuration options related to the man page output
-#---------------------------------------------------------------------------
-
-# If the GENERATE_MAN tag is set to YES, doxygen will generate man pages for
-# classes and files.
-# The default value is: NO.
-
-GENERATE_MAN           = NO
-
-# The MAN_OUTPUT tag is used to specify where the man pages will be put. If a
-# relative path is entered the value of OUTPUT_DIRECTORY will be put in front of
-# it. A directory man3 will be created inside the directory specified by
-# MAN_OUTPUT.
-# The default directory is: man.
-# This tag requires that the tag GENERATE_MAN is set to YES.
-
-MAN_OUTPUT             = man
-
-# The MAN_EXTENSION tag determines the extension that is added to the generated
-# man pages. In case the manual section does not start with a number, the number
-# 3 is prepended. The dot (.) at the beginning of the MAN_EXTENSION tag is
-# optional.
-# The default value is: .3.
-# This tag requires that the tag GENERATE_MAN is set to YES.
-
-MAN_EXTENSION          = .3
-
-# The MAN_SUBDIR tag determines the name of the directory created within
-# MAN_OUTPUT in which the man pages are placed. If defaults to man followed by
-# MAN_EXTENSION with the initial . removed.
-# This tag requires that the tag GENERATE_MAN is set to YES.
-
-MAN_SUBDIR             =
-
-# If the MAN_LINKS tag is set to YES and doxygen generates man output, then it
-# will generate one additional man file for each entity documented in the real
-# man page(s). These additional files only source the real man page, but without
-# them the man command would be unable to find the correct page.
-# The default value is: NO.
-# This tag requires that the tag GENERATE_MAN is set to YES.
-
-MAN_LINKS              = NO
-
-#---------------------------------------------------------------------------
-# Configuration options related to the XML output
-#---------------------------------------------------------------------------
-
-# If the GENERATE_XML tag is set to YES, doxygen will generate an XML file that
-# captures the structure of the code including all documentation.
-# The default value is: NO.
-
-GENERATE_XML           = YES
-
-# The XML_OUTPUT tag is used to specify where the XML pages will be put. If a
-# relative path is entered the value of OUTPUT_DIRECTORY will be put in front of
-# it.
-# The default directory is: xml.
-# This tag requires that the tag GENERATE_XML is set to YES.
-
-XML_OUTPUT             = xml
-
-# If the XML_PROGRAMLISTING tag is set to YES, doxygen will dump the program
-# listings (including syntax highlighting and cross-referencing information) to
-# the XML output. Note that enabling this will significantly increase the size
-# of the XML output.
-# The default value is: YES.
-# This tag requires that the tag GENERATE_XML is set to YES.
-
-XML_PROGRAMLISTING     = YES
-
-#---------------------------------------------------------------------------
-# Configuration options related to the DOCBOOK output
-#---------------------------------------------------------------------------
-
-# If the GENERATE_DOCBOOK tag is set to YES, doxygen will generate Docbook files
-# that can be used to generate PDF.
-# The default value is: NO.
-
-GENERATE_DOCBOOK       = NO
-
-# The DOCBOOK_OUTPUT tag is used to specify where the Docbook pages will be put.
-# If a relative path is entered the value of OUTPUT_DIRECTORY will be put in
-# front of it.
-# The default directory is: docbook.
-# This tag requires that the tag GENERATE_DOCBOOK is set to YES.
-
-DOCBOOK_OUTPUT         = docbook
-
-# If the DOCBOOK_PROGRAMLISTING tag is set to YES, doxygen will include the
-# program listings (including syntax highlighting and cross-referencing
-# information) to the DOCBOOK output. Note that enabling this will significantly
-# increase the size of the DOCBOOK output.
-# The default value is: NO.
-# This tag requires that the tag GENERATE_DOCBOOK is set to YES.
-
-DOCBOOK_PROGRAMLISTING = NO
-
-#---------------------------------------------------------------------------
-# Configuration options for the AutoGen Definitions output
-#---------------------------------------------------------------------------
-
-# If the GENERATE_AUTOGEN_DEF tag is set to YES, doxygen will generate an
-# AutoGen Definitions (see http://autogen.sourceforge.net/) file that captures
-# the structure of the code including all documentation. Note that this feature
-# is still experimental and incomplete at the moment.
-# The default value is: NO.
-
-GENERATE_AUTOGEN_DEF   = NO
-
-#---------------------------------------------------------------------------
-# Configuration options related to the Perl module output
-#---------------------------------------------------------------------------
-
-# If the GENERATE_PERLMOD tag is set to YES, doxygen will generate a Perl module
-# file that captures the structure of the code including all documentation.
-#
-# Note that this feature is still experimental and incomplete at the moment.
-# The default value is: NO.
-
-GENERATE_PERLMOD       = NO
-
-# If the PERLMOD_LATEX tag is set to YES, doxygen will generate the necessary
-# Makefile rules, Perl scripts and LaTeX code to be able to generate PDF and DVI
-# output from the Perl module output.
-# The default value is: NO.
-# This tag requires that the tag GENERATE_PERLMOD is set to YES.
-
-PERLMOD_LATEX          = NO
-
-# If the PERLMOD_PRETTY tag is set to YES, the Perl module output will be nicely
-# formatted so it can be parsed by a human reader. This is useful if you want to
-# understand what is going on. On the other hand, if this tag is set to NO, the
-# size of the Perl module output will be much smaller and Perl will parse it
-# just the same.
-# The default value is: YES.
-# This tag requires that the tag GENERATE_PERLMOD is set to YES.
-
-PERLMOD_PRETTY         = YES
-
-# The names of the make variables in the generated doxyrules.make file are
-# prefixed with the string contained in PERLMOD_MAKEVAR_PREFIX. This is useful
-# so different doxyrules.make files included by the same Makefile don't
-# overwrite each other's variables.
-# This tag requires that the tag GENERATE_PERLMOD is set to YES.
-
-PERLMOD_MAKEVAR_PREFIX =
 
 #---------------------------------------------------------------------------
 # Configuration options related to the preprocessor
@@ -2162,321 +1889,29 @@ EXTERNAL_PAGES         = YES
 
 PERL_PATH              = /usr/bin/perl
 
+
 #---------------------------------------------------------------------------
-# Configuration options related to the dot tool
+# Configuration options related to the XML output
 #---------------------------------------------------------------------------
-
-# If the CLASS_DIAGRAMS tag is set to YES, doxygen will generate a class diagram
-# (in HTML and LaTeX) for classes with base or super classes. Setting the tag to
-# NO turns the diagrams off. Note that this option also works with HAVE_DOT
-# disabled, but it is recommended to install and use dot, since it yields more
-# powerful graphs.
-# The default value is: YES.
-
-CLASS_DIAGRAMS         = YES
-
-# You can define message sequence charts within doxygen comments using the \msc
-# command. Doxygen will then run the mscgen tool (see:
-# http://www.mcternan.me.uk/mscgen/)) to produce the chart and insert it in the
-# documentation. The MSCGEN_PATH tag allows you to specify the directory where
-# the mscgen tool resides. If left empty the tool is assumed to be found in the
-# default search path.
-
-MSCGEN_PATH            =
-
-# You can include diagrams made with dia in doxygen documentation. Doxygen will
-# then run dia to produce the diagram and insert it in the documentation. The
-# DIA_PATH tag allows you to specify the directory where the dia binary resides.
-# If left empty dia is assumed to be found in the default search path.
-
-DIA_PATH               =
-
-# If set to YES the inheritance and collaboration graphs will hide inheritance
-# and usage relations if the target is undocumented or is not a class.
-# The default value is: YES.
-
-HIDE_UNDOC_RELATIONS   = YES
-
-# If you set the HAVE_DOT tag to YES then doxygen will assume the dot tool is
-# available from the path. This tool is part of Graphviz (see:
-# http://www.graphviz.org/), a graph visualization toolkit from AT&T and Lucent
-# Bell Labs. The other options in this section have no effect if this option is
-# set to NO
+# If the GENERATE_XML tag is set to YES, doxygen will generate an XML file that
+# captures the structure of the code including all documentation.
 # The default value is: NO.
 
-HAVE_DOT               = NO
+GENERATE_XML           = YES
 
-# The DOT_NUM_THREADS specifies the number of dot invocations doxygen is allowed
-# to run in parallel. When set to 0 doxygen will base this on the number of
-# processors available in the system. You can set it explicitly to a value
-# larger than 0 to get control over the balance between CPU load and processing
-# speed.
-# Minimum value: 0, maximum value: 32, default value: 0.
-# This tag requires that the tag HAVE_DOT is set to YES.
+# The XML_OUTPUT tag is used to specify where the XML pages will be put. If a
+# relative path is entered the value of OUTPUT_DIRECTORY will be put in front of
+# it.
+# The default directory is: xml.
+# This tag requires that the tag GENERATE_XML is set to YES.
 
-DOT_NUM_THREADS        = 0
+XML_OUTPUT             = xml
 
-# When you want a differently looking font in the dot files that doxygen
-# generates you can specify the font name using DOT_FONTNAME. You need to make
-# sure dot is able to find the font, which can be done by putting it in a
-# standard location or by setting the DOTFONTPATH environment variable or by
-# setting DOT_FONTPATH to the directory containing the font.
-# The default value is: Helvetica.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-DOT_FONTNAME           = Helvetica
-
-# The DOT_FONTSIZE tag can be used to set the size (in points) of the font of
-# dot graphs.
-# Minimum value: 4, maximum value: 24, default value: 10.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-DOT_FONTSIZE           = 10
-
-# By default doxygen will tell dot to use the default font as specified with
-# DOT_FONTNAME. If you specify a different font using DOT_FONTNAME you can set
-# the path where dot can find it using this tag.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-DOT_FONTPATH           =
-
-# If the CLASS_GRAPH tag is set to YES then doxygen will generate a graph for
-# each documented class showing the direct and indirect inheritance relations.
-# Setting this tag to YES will force the CLASS_DIAGRAMS tag to NO.
+# If the XML_PROGRAMLISTING tag is set to YES, doxygen will dump the program
+# listings (including syntax highlighting and cross-referencing information) to
+# the XML output. Note that enabling this will significantly increase the size
+# of the XML output.
 # The default value is: YES.
-# This tag requires that the tag HAVE_DOT is set to YES.
+# This tag requires that the tag GENERATE_XML is set to YES.
 
-CLASS_GRAPH            = YES
-
-# If the COLLABORATION_GRAPH tag is set to YES then doxygen will generate a
-# graph for each documented class showing the direct and indirect implementation
-# dependencies (inheritance, containment, and class references variables) of the
-# class with other documented classes.
-# The default value is: YES.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-COLLABORATION_GRAPH    = YES
-
-# If the GROUP_GRAPHS tag is set to YES then doxygen will generate a graph for
-# groups, showing the direct groups dependencies.
-# The default value is: YES.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-GROUP_GRAPHS           = YES
-
-# If the UML_LOOK tag is set to YES, doxygen will generate inheritance and
-# collaboration diagrams in a style similar to the OMG's Unified Modeling
-# Language.
-# The default value is: NO.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-UML_LOOK               = NO
-
-# If the UML_LOOK tag is enabled, the fields and methods are shown inside the
-# class node. If there are many fields or methods and many nodes the graph may
-# become too big to be useful. The UML_LIMIT_NUM_FIELDS threshold limits the
-# number of items for each type to make the size more manageable. Set this to 0
-# for no limit. Note that the threshold may be exceeded by 50% before the limit
-# is enforced. So when you set the threshold to 10, up to 15 fields may appear,
-# but if the number exceeds 15, the total amount of fields shown is limited to
-# 10.
-# Minimum value: 0, maximum value: 100, default value: 10.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-UML_LIMIT_NUM_FIELDS   = 10
-
-# If the TEMPLATE_RELATIONS tag is set to YES then the inheritance and
-# collaboration graphs will show the relations between templates and their
-# instances.
-# The default value is: NO.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-TEMPLATE_RELATIONS     = NO
-
-# If the INCLUDE_GRAPH, ENABLE_PREPROCESSING and SEARCH_INCLUDES tags are set to
-# YES then doxygen will generate a graph for each documented file showing the
-# direct and indirect include dependencies of the file with other documented
-# files.
-# The default value is: YES.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-INCLUDE_GRAPH          = YES
-
-# If the INCLUDED_BY_GRAPH, ENABLE_PREPROCESSING and SEARCH_INCLUDES tags are
-# set to YES then doxygen will generate a graph for each documented file showing
-# the direct and indirect include dependencies of the file with other documented
-# files.
-# The default value is: YES.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-INCLUDED_BY_GRAPH      = YES
-
-# If the CALL_GRAPH tag is set to YES then doxygen will generate a call
-# dependency graph for every global function or class method.
-#
-# Note that enabling this option will significantly increase the time of a run.
-# So in most cases it will be better to enable call graphs for selected
-# functions only using the \callgraph command. Disabling a call graph can be
-# accomplished by means of the command \hidecallgraph.
-# The default value is: NO.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-CALL_GRAPH             = NO
-
-# If the CALLER_GRAPH tag is set to YES then doxygen will generate a caller
-# dependency graph for every global function or class method.
-#
-# Note that enabling this option will significantly increase the time of a run.
-# So in most cases it will be better to enable caller graphs for selected
-# functions only using the \callergraph command. Disabling a caller graph can be
-# accomplished by means of the command \hidecallergraph.
-# The default value is: NO.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-CALLER_GRAPH           = NO
-
-# If the GRAPHICAL_HIERARCHY tag is set to YES then doxygen will graphical
-# hierarchy of all classes instead of a textual one.
-# The default value is: YES.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-GRAPHICAL_HIERARCHY    = YES
-
-# If the DIRECTORY_GRAPH tag is set to YES then doxygen will show the
-# dependencies a directory has on other directories in a graphical way. The
-# dependency relations are determined by the #include relations between the
-# files in the directories.
-# The default value is: YES.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-DIRECTORY_GRAPH        = YES
-
-# The DOT_IMAGE_FORMAT tag can be used to set the image format of the images
-# generated by dot. For an explanation of the image formats see the section
-# output formats in the documentation of the dot tool (Graphviz (see:
-# http://www.graphviz.org/)).
-# Note: If you choose svg you need to set HTML_FILE_EXTENSION to xhtml in order
-# to make the SVG files visible in IE 9+ (other browsers do not have this
-# requirement).
-# Possible values are: png, jpg, gif, svg, png:gd, png:gd:gd, png:cairo,
-# png:cairo:gd, png:cairo:cairo, png:cairo:gdiplus, png:gdiplus and
-# png:gdiplus:gdiplus.
-# The default value is: png.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-DOT_IMAGE_FORMAT       = png
-
-# If DOT_IMAGE_FORMAT is set to svg, then this option can be set to YES to
-# enable generation of interactive SVG images that allow zooming and panning.
-#
-# Note that this requires a modern browser other than Internet Explorer. Tested
-# and working are Firefox, Chrome, Safari, and Opera.
-# Note: For IE 9+ you need to set HTML_FILE_EXTENSION to xhtml in order to make
-# the SVG files visible. Older versions of IE do not have SVG support.
-# The default value is: NO.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-INTERACTIVE_SVG        = NO
-
-# The DOT_PATH tag can be used to specify the path where the dot tool can be
-# found. If left blank, it is assumed the dot tool can be found in the path.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-DOT_PATH               =
-
-# The DOTFILE_DIRS tag can be used to specify one or more directories that
-# contain dot files that are included in the documentation (see the \dotfile
-# command).
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-DOTFILE_DIRS           =
-
-# The MSCFILE_DIRS tag can be used to specify one or more directories that
-# contain msc files that are included in the documentation (see the \mscfile
-# command).
-
-MSCFILE_DIRS           =
-
-# The DIAFILE_DIRS tag can be used to specify one or more directories that
-# contain dia files that are included in the documentation (see the \diafile
-# command).
-
-DIAFILE_DIRS           =
-
-# When using plantuml, the PLANTUML_JAR_PATH tag should be used to specify the
-# path where java can find the plantuml.jar file. If left blank, it is assumed
-# PlantUML is not used or called during a preprocessing step. Doxygen will
-# generate a warning when it encounters a \startuml command in this case and
-# will not generate output for the diagram.
-
-PLANTUML_JAR_PATH      =
-
-# When using plantuml, the PLANTUML_CFG_FILE tag can be used to specify a
-# configuration file for plantuml.
-
-PLANTUML_CFG_FILE      =
-
-# When using plantuml, the specified paths are searched for files specified by
-# the !include statement in a plantuml block.
-
-PLANTUML_INCLUDE_PATH  =
-
-# The DOT_GRAPH_MAX_NODES tag can be used to set the maximum number of nodes
-# that will be shown in the graph. If the number of nodes in a graph becomes
-# larger than this value, doxygen will truncate the graph, which is visualized
-# by representing a node as a red box. Note that doxygen if the number of direct
-# children of the root node in a graph is already larger than
-# DOT_GRAPH_MAX_NODES then the graph will not be shown at all. Also note that
-# the size of a graph can be further restricted by MAX_DOT_GRAPH_DEPTH.
-# Minimum value: 0, maximum value: 10000, default value: 50.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-DOT_GRAPH_MAX_NODES    = 50
-
-# The MAX_DOT_GRAPH_DEPTH tag can be used to set the maximum depth of the graphs
-# generated by dot. A depth value of 3 means that only nodes reachable from the
-# root by following a path via at most 3 edges will be shown. Nodes that lay
-# further from the root node will be omitted. Note that setting this option to 1
-# or 2 may greatly reduce the computation time needed for large code bases. Also
-# note that the size of a graph can be further restricted by
-# DOT_GRAPH_MAX_NODES. Using a depth of 0 means no depth restriction.
-# Minimum value: 0, maximum value: 1000, default value: 0.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-MAX_DOT_GRAPH_DEPTH    = 0
-
-# Set the DOT_TRANSPARENT tag to YES to generate images with a transparent
-# background. This is disabled by default, because dot on Windows does not seem
-# to support this out of the box.
-#
-# Warning: Depending on the platform used, enabling this option may lead to
-# badly anti-aliased labels on the edges of a graph (i.e. they become hard to
-# read).
-# The default value is: NO.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-DOT_TRANSPARENT        = NO
-
-# Set the DOT_MULTI_TARGETS tag to YES to allow dot to generate multiple output
-# files in one run (i.e. multiple -o and -T options on the command line). This
-# makes dot run faster, but since only newer versions of dot (>1.8.10) support
-# this, this feature is disabled by default.
-# The default value is: NO.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-DOT_MULTI_TARGETS      = NO
-
-# If the GENERATE_LEGEND tag is set to YES doxygen will generate a legend page
-# explaining the meaning of the various boxes and arrows in the dot generated
-# graphs.
-# The default value is: YES.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-GENERATE_LEGEND        = YES
-
-# If the DOT_CLEANUP tag is set to YES, doxygen will remove the intermediate dot
-# files that are used to generate the various graphs.
-# The default value is: YES.
-# This tag requires that the tag HAVE_DOT is set to YES.
-
-DOT_CLEANUP            = YES
+XML_PROGRAMLISTING     = YES

docs/conf.py.in

@@ -12,10 +12,8 @@
 #
 import os
 import sys
-# sys.path.insert(0, os.path.abspath('.'))
-sys.path.insert(0, os.path.abspath('../bin/'))
-#sys.path.insert(0, '/home/l_frojdh/sls/build/bin')
-#sys.path.insert(0, @CMAKE_CURRENT_BINARY_DIR@)
+sys.path.insert(0, os.path.abspath('..'))
+
 print(sys.path)
 
 # -- Project information -----------------------------------------------------
@@ -31,7 +29,6 @@ version = '@PROJECT_VERSION@'
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
 # ones.
 extensions = ['breathe',
-              'sphinx_rtd_theme',
               'sphinx.ext.autodoc',
               'sphinx.ext.napoleon',
              ]

docs/src/ClusterFile.rst

@@ -0,0 +1,7 @@
+ClusterFile
+=============
+
+.. doxygenclass:: aare::ClusterFile
+   :members:
+   :undoc-members:
+   :private-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/File.rst

@@ -4,4 +4,5 @@ File
 
 .. doxygenclass:: aare::File
    :members:
-   :undoc-members:
+   :undoc-members:
+   :private-members:

docs/src/Frame.rst

@@ -4,4 +4,5 @@ Frame
 
 .. doxygenclass:: aare::Frame
    :members:
-   :undoc-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/JungfrauDataFile.rst

@@ -0,0 +1,25 @@
+JungfrauDataFile
+==================
+
+JungfrauDataFile is a class to read the .dat files that are produced by Aldo's receiver. 
+It is mostly used for calibration. 
+
+The structure of the file is:
+
+* JungfrauDataHeader
+* Binary data (256x256, 256x1024 or 512x1024)
+* JungfrauDataHeader
+* ...
+
+There is no metadata indicating number of frames or the size of the image, but this
+will be infered by this reader. 
+
+.. doxygenstruct:: aare::JungfrauDataHeader
+   :members:
+   :undoc-members:
+   :private-members:
+
+.. doxygenclass:: aare::JungfrauDataFile
+   :members:
+   :undoc-members:
+   :private-members:

docs/src/Pedestal.rst

@@ -4,4 +4,5 @@ Pedestal
 
 .. doxygenclass:: aare::Pedestal
    :members:
-   :undoc-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/Tests.rst

@@ -0,0 +1,47 @@
+****************
+Tests
+****************
+
+We test the code both from the C++ and Python API. By default only tests that does not require image data is run. 
+
+C++
+~~~~~~~~~~~~~~~~~~
+
+.. code-block:: bash
+
+    mkdir build
+    cd build
+    cmake .. -DAARE_TESTS=ON
+    make -j 4
+
+    export AARE_TEST_DATA=/path/to/test/data
+    ./run_test [.files] #or using ctest, [.files] is the option to include tests needing data
+
+
+
+Python
+~~~~~~~~~~~~~~~~~~
+
+.. code-block:: bash
+
+    #From the root dir of the library
+    python -m pytest python/tests --files # passing --files will run the tests needing data
+
+
+
+Getting the test data
+~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. attention ::
+
+    The tests needing the test data are not run by default. To make the data available, you need to set the environment variable
+    AARE_TEST_DATA to the path of the test data directory. Then pass either [.files] for the C++ tests or --files for Python
+
+The image files needed for the test are large and are not included in the repository. They are stored
+using GIT LFS in a separate repository. To get the test data, you need to clone the repository.
+To do this, you need to have GIT LFS installed. You can find instructions on how to install it here: https://git-lfs.github.com/
+Once you have GIT LFS installed, you can clone the repository like any normal repo using:
+
+.. code-block:: bash
+    
+    git clone https://gitea.psi.ch/detectors/aare-test-data.git

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/algorithm.rst

@@ -0,0 +1,5 @@
+algorithm
+=============
+
+.. doxygenfile:: algorithm.hpp
+

docs/src/index.rst

@@ -2,24 +2,65 @@ 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>`_
+
+
 
-    Hello
 
 .. toctree::
-    :caption: C++ API
-    :maxdepth: 1
+    :caption: Installation
+    :maxdepth: 3
     
-    NDArray
-    NDView
-    Frame
-    File
-    Dtype
-    ClusterFinder
-    Pedestal
-    VarClusterFinder
+    Installation
+    Requirements
+    Consume
 
 .. toctree::
     :caption: Python API
     :maxdepth: 1
     
     pyFile
+    pyCtbRawFile
+    pyClusterFile
+    pyClusterVector
+    pyJungfrauDataFile
+    pyRawFile
+    pyRawMasterFile
+    pyVarClusterFinder
+
+    pyFit
+
+
+.. toctree::
+    :caption: C++ API
+    :maxdepth: 1
+    
+    algorithm
+    NDArray
+    NDView
+    Frame
+    File
+    Dtype
+    ClusterFinder
+    ClusterFinderMT
+    ClusterFile
+    ClusterVector
+    JungfrauDataFile
+    Pedestal
+    RawFile
+    RawSubFile
+    RawMasterFile
+    VarClusterFinder
+
+
+
+.. toctree::
+    :caption: Developer
+    :maxdepth: 3
+    
+    Tests

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/pyFit.rst

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

docs/src/pyJungfrauDataFile.rst

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

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:

etc/dev-env.yml

@@ -0,0 +1,15 @@
+name: dev-environment
+channels:
+  - conda-forge
+dependencies:
+  - anaconda-client
+  - doxygen 
+  - sphinx=7.1.2 
+  - breathe 
+  - pybind11 
+  - sphinx_rtd_theme 
+  - furo 
+  - nlohmann_json 
+  - zeromq 
+  - fmt 
+  - numpy

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/CalculateEta.hpp

@@ -0,0 +1,122 @@
+#pragma once
+
+#include "aare/Cluster.hpp"
+#include "aare/ClusterVector.hpp"
+#include "aare/NDArray.hpp"
+
+namespace aare {
+
+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;
+
+template <typename T> struct Eta2 {
+    double x;
+    double y;
+    int c;
+    T sum;
+};
+
+/**
+ * @brief Calculate the eta2 values for all clusters in a Clsutervector
+ */
+template <typename ClusterType,
+          typename = std::enable_if_t<is_cluster_v<ClusterType>>>
+NDArray<double, 2> calculate_eta2(const ClusterVector<ClusterType> &clusters) {
+    NDArray<double, 2> eta2({static_cast<int64_t>(clusters.size()), 2});
+
+    for (size_t i = 0; i < clusters.size(); i++) {
+        auto e = calculate_eta2(clusters.at(i));
+        eta2(i, 0) = e.x;
+        eta2(i, 1) = e.y;
+    }
+
+    return eta2;
+}
+
+/**
+ * @brief Calculate the eta2 values for a generic sized cluster and return them
+ * in a Eta2 struct containing etay, etax and the index of the respective 2x2
+ * subcluster.
+ */
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType>
+Eta2<T>
+calculate_eta2(const Cluster<T, ClusterSizeX, ClusterSizeY, CoordType> &cl) {
+    Eta2<T> eta{};
+
+    auto max_sum = cl.max_sum_2x2();
+    eta.sum = max_sum.first;
+    auto c = max_sum.second;
+
+    size_t index_bottom_left_max_2x2_subcluster =
+        (int(c / (ClusterSizeX - 1))) * ClusterSizeX + c % (ClusterSizeX - 1);
+
+    if ((cl.data[index_bottom_left_max_2x2_subcluster] +
+         cl.data[index_bottom_left_max_2x2_subcluster + 1]) != 0)
+        eta.x = static_cast<double>(
+                    cl.data[index_bottom_left_max_2x2_subcluster + 1]) /
+                static_cast<double>(
+                    (cl.data[index_bottom_left_max_2x2_subcluster] +
+                     cl.data[index_bottom_left_max_2x2_subcluster + 1]));
+
+    if ((cl.data[index_bottom_left_max_2x2_subcluster] +
+         cl.data[index_bottom_left_max_2x2_subcluster + ClusterSizeX]) != 0)
+        eta.y =
+            static_cast<double>(
+                cl.data[index_bottom_left_max_2x2_subcluster + ClusterSizeX]) /
+            static_cast<double>(
+                (cl.data[index_bottom_left_max_2x2_subcluster] +
+                 cl.data[index_bottom_left_max_2x2_subcluster + ClusterSizeX]));
+
+    eta.c = c; // TODO only supported for 2x2 and 3x3 clusters -> at least no
+               // underyling enum class
+    return eta;
+}
+
+// calculates Eta3 for 3x3 cluster based on code from analyze_cluster
+// TODO only supported for 3x3 Clusters
+template <typename T> Eta2<T> calculate_eta3(const Cluster<T, 3, 3> &cl) {
+
+    Eta2<T> eta{};
+
+    T sum = 0;
+
+    std::for_each(std::begin(cl.data), std::end(cl.data),
+                  [&sum](T x) { sum += x; });
+
+    eta.sum = sum;
+
+    eta.c = corner::cBottomLeft;
+
+    if ((cl.data[3] + cl.data[4] + cl.data[5]) != 0)
+
+        eta.x = static_cast<double>(-cl.data[3] + cl.data[3 + 2]) /
+
+                (cl.data[3] + cl.data[4] + cl.data[5]);
+
+    if ((cl.data[1] + cl.data[4] + cl.data[7]) != 0)
+
+        eta.y = static_cast<double>(-cl.data[1] + cl.data[2 * 3 + 1]) /
+
+                (cl.data[1] + cl.data[4] + cl.data[7]);
+
+    return eta;
+}
+
+} // 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/Cluster.hpp

@@ -0,0 +1,121 @@
+
+/************************************************
+ * @file Cluster.hpp
+ * @short definition of cluster, where CoordType (x,y) give
+ * the cluster center coordinates and data the actual cluster data
+ * cluster size is given as template parameters
+ ***********************************************/
+
+#pragma once
+
+#include <algorithm>
+#include <array>
+#include <cstdint>
+#include <numeric>
+#include <type_traits>
+
+namespace aare {
+
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = int16_t>
+constexpr bool is_valid_cluster =
+    std::is_arithmetic_v<T> && std::is_integral_v<CoordType> &&
+    (ClusterSizeX > 0) && (ClusterSizeY > 0);
+
+// requires clause c++20 maybe update
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType = int16_t,
+          typename Enable = std::enable_if_t<
+              is_valid_cluster<T, ClusterSizeX, ClusterSizeY, CoordType>>>
+struct Cluster {
+    CoordType x;
+    CoordType y;
+    T data[ClusterSizeX * ClusterSizeY];
+
+    T sum() const {
+        return std::accumulate(data, data + ClusterSizeX * ClusterSizeY, 0);
+    }
+
+    std::pair<T, int> max_sum_2x2() const {
+
+        constexpr size_t num_2x2_subclusters =
+            (ClusterSizeX - 1) * (ClusterSizeY - 1);
+
+        std::array<T, num_2x2_subclusters> sum_2x2_subcluster;
+        for (size_t i = 0; i < ClusterSizeY - 1; ++i) {
+            for (size_t j = 0; j < ClusterSizeX - 1; ++j)
+                sum_2x2_subcluster[i * (ClusterSizeX - 1) + j] =
+                    data[i * ClusterSizeX + j] +
+                    data[i * ClusterSizeX + j + 1] +
+                    data[(i + 1) * ClusterSizeX + j] +
+                    data[(i + 1) * ClusterSizeX + j + 1];
+        }
+
+        int index = std::max_element(sum_2x2_subcluster.begin(),
+                                     sum_2x2_subcluster.end()) -
+                    sum_2x2_subcluster.begin();
+        return std::make_pair(sum_2x2_subcluster[index], index);
+    }
+};
+
+// Specialization for 2x2 clusters (only one sum exists)
+template <typename T> struct Cluster<T, 2, 2, int16_t> {
+    int16_t x;
+    int16_t y;
+    T data[4];
+
+    T sum() const { return std::accumulate(data, data + 4, 0); }
+
+    std::pair<T, int> max_sum_2x2() const {
+        return std::make_pair(data[0] + data[1] + data[2] + data[3],
+                              0); // Only one possible 2x2 sum
+    }
+};
+
+// Specialization for 3x3 clusters
+template <typename T> struct Cluster<T, 3, 3, int16_t> {
+    int16_t x;
+    int16_t y;
+    T data[9];
+
+    T sum() const { return std::accumulate(data, data + 9, 0); }
+
+    std::pair<T, int> max_sum_2x2() const {
+        std::array<T, 4> sum_2x2_subclusters;
+        sum_2x2_subclusters[0] = data[0] + data[1] + data[3] + data[4];
+        sum_2x2_subclusters[1] = data[1] + data[2] + data[4] + data[5];
+        sum_2x2_subclusters[2] = data[3] + data[4] + data[6] + data[7];
+        sum_2x2_subclusters[3] = data[4] + data[5] + data[7] + data[8];
+        int index = std::max_element(sum_2x2_subclusters.begin(),
+                                     sum_2x2_subclusters.end()) -
+                    sum_2x2_subclusters.begin();
+        return std::make_pair(sum_2x2_subclusters[index], index);
+    }
+};
+
+// Type Traits for is_cluster_type
+template <typename T>
+struct is_cluster : std::false_type {}; // Default case: Not a Cluster
+
+template <typename T, uint8_t X, uint8_t Y, typename CoordType>
+struct is_cluster<Cluster<T, X, Y, CoordType>> : std::true_type {}; // Cluster
+
+template <typename T> constexpr bool is_cluster_v = is_cluster<T>::value;
+
+template <typename ClusterType,
+          typename = std::enable_if_t<is_cluster_v<ClusterType>>>
+struct extract_template_arguments; // Forward declaration
+
+// helper struct to extract template argument
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType>
+struct extract_template_arguments<
+    Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
+
+    using value_type = T;
+    static constexpr int cluster_size_x = ClusterSizeX;
+    static constexpr int cluster_size_y = ClusterSizeY;
+    using coordtype = CoordType;
+};
+
+} // namespace aare

include/aare/ClusterCollector.hpp

@@ -0,0 +1,54 @@
+#pragma once
+#include <atomic>
+#include <thread>
+
+#include "aare/ClusterFinderMT.hpp"
+#include "aare/ClusterVector.hpp"
+#include "aare/ProducerConsumerQueue.hpp"
+
+namespace aare {
+
+template <typename ClusterType,
+          typename = std::enable_if_t<is_cluster_v<ClusterType>>>
+class ClusterCollector {
+    ProducerConsumerQueue<ClusterVector<ClusterType>> *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<ClusterType>> m_clusters;
+
+    void process() {
+        m_stopped = false;
+        fmt::print("ClusterCollector started\n");
+        while (!m_stop_requested || !m_source->isEmpty()) {
+            if (ClusterVector<ClusterType> *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<ClusterType, uint16_t, double> *source) {
+        m_source = source->sink();
+        m_thread = std::thread(&ClusterCollector::process, this);
+    }
+    void stop() {
+        m_stop_requested = true;
+        m_thread.join();
+    }
+    std::vector<ClusterVector<ClusterType>> 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,454 @@
+#pragma once
+
+#include "aare/Cluster.hpp"
+#include "aare/ClusterVector.hpp"
+#include "aare/GainMap.hpp"
+#include "aare/NDArray.hpp"
+#include "aare/defs.hpp"
+#include <filesystem>
+#include <fstream>
+#include <optional>
+
+namespace aare {
+
+/*
+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
+....
+*/
+
+// TODO: change to support any type of clusters, e.g. header line with
+// clsuter_size_x, cluster_size_y,
+/**
+ * @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] * number_of_clusters
+ *       int32_t frame_number
+ *       uint32_t number_of_clusters
+ *       etc.
+ */
+template <typename ClusterType,
+          typename Enable = std::enable_if_t<is_cluster_v<ClusterType>>>
+class ClusterFile {
+    FILE *fp{};
+    uint32_t m_num_left{};    /*Number of photons left in frame*/
+    size_t m_chunk_size{};    /*Number of clusters to read at a time*/
+    const std::string m_mode; /*Mode to open the file in*/
+    std::optional<ROI> m_roi; /*Region of interest, will be applied if set*/
+    std::optional<NDArray<int32_t, 2>>
+        m_noise_map; /*Noise map to cut photons, will be applied if set*/
+    std::optional<GainMap> m_gain_map; /*Gain map to apply to the clusters, will
+                                          be applied if set*/
+
+  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<ClusterType> 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<ClusterType> read_frame();
+
+    void write_frame(const ClusterVector<ClusterType> &clusters);
+
+    /**
+     * @brief Return the chunk size
+     */
+    size_t chunk_size() const { return m_chunk_size; }
+
+    /**
+     * @brief Set the region of interest to use when reading clusters. If set
+     * only clusters within the ROI will be read.
+     */
+    void set_roi(ROI roi);
+
+    /**
+     * @brief Set the noise map to use when reading clusters. If set clusters
+     * below the noise level will be discarded. Selection criteria one of:
+     * Central pixel above noise, highest 2x2 sum above 2 * noise, total sum
+     * above 3 * noise.
+     */
+    void set_noise_map(const NDView<int32_t, 2> noise_map);
+
+    /**
+     * @brief Set the gain map to use when reading clusters. If set the gain map
+     * will be applied to the clusters that pass ROI and noise_map selection.
+     */
+    void set_gain_map(const NDView<double, 2> gain_map);
+
+    void set_gain_map(const GainMap &gain_map);
+
+    void set_gain_map(const GainMap &&gain_map);
+
+    /**
+     * @brief Close the file. If not closed the file will be closed in the
+     * destructor
+     */
+    void close();
+
+  private:
+    ClusterVector<ClusterType> read_clusters_with_cut(size_t n_clusters);
+    ClusterVector<ClusterType> read_clusters_without_cut(size_t n_clusters);
+    ClusterVector<ClusterType> read_frame_with_cut();
+    ClusterVector<ClusterType> read_frame_without_cut();
+    bool is_selected(ClusterType &cl);
+    ClusterType read_one_cluster();
+};
+
+template <typename ClusterType, typename Enable>
+ClusterFile<ClusterType, Enable>::ClusterFile(
+    const std::filesystem::path &fname, size_t chunk_size,
+    const std::string &mode)
+    : m_chunk_size(chunk_size), m_mode(mode) {
+
+    if (mode == "r") {
+        fp = fopen(fname.c_str(), "rb");
+        if (!fp) {
+            throw std::runtime_error("Could not open file for reading: " +
+                                     fname.string());
+        }
+    } else if (mode == "w") {
+        fp = fopen(fname.c_str(), "wb");
+        if (!fp) {
+            throw std::runtime_error("Could not open file for writing: " +
+                                     fname.string());
+        }
+    } else if (mode == "a") {
+        fp = fopen(fname.c_str(), "ab");
+        if (!fp) {
+            throw std::runtime_error("Could not open file for appending: " +
+                                     fname.string());
+        }
+    } else {
+        throw std::runtime_error("Unsupported mode: " + mode);
+    }
+}
+
+template <typename ClusterType, typename Enable>
+ClusterFile<ClusterType, Enable>::~ClusterFile() {
+    close();
+}
+
+template <typename ClusterType, typename Enable>
+void ClusterFile<ClusterType, Enable>::close() {
+    if (fp) {
+        fclose(fp);
+        fp = nullptr;
+    }
+}
+template <typename ClusterType, typename Enable>
+void ClusterFile<ClusterType, Enable>::set_roi(ROI roi) {
+    m_roi = roi;
+}
+template <typename ClusterType, typename Enable>
+void ClusterFile<ClusterType, Enable>::set_noise_map(
+    const NDView<int32_t, 2> noise_map) {
+    m_noise_map = NDArray<int32_t, 2>(noise_map);
+}
+template <typename ClusterType, typename Enable>
+void ClusterFile<ClusterType, Enable>::set_gain_map(
+    const NDView<double, 2> gain_map) {
+    m_gain_map = GainMap(gain_map);
+}
+
+template <typename ClusterType, typename Enable>
+void ClusterFile<ClusterType, Enable>::set_gain_map(const GainMap &gain_map) {
+    m_gain_map = gain_map;
+}
+
+template <typename ClusterType, typename Enable>
+void ClusterFile<ClusterType, Enable>::set_gain_map(const GainMap &&gain_map) {
+    m_gain_map = gain_map;
+}
+
+// TODO generally supported for all clsuter types
+template <typename ClusterType, typename Enable>
+void ClusterFile<ClusterType, Enable>::write_frame(
+    const ClusterVector<ClusterType> &clusters) {
+    if (m_mode != "w" && m_mode != "a") {
+        throw std::runtime_error("File not opened for writing");
+    }
+    if (!(clusters.cluster_size_x() == 3) &&
+        !(clusters.cluster_size_y() == 3)) {
+        throw std::runtime_error("Only 3x3 clusters are supported");
+    }
+    int32_t frame_number = clusters.frame_number();
+    fwrite(&frame_number, sizeof(frame_number), 1, fp);
+    uint32_t n_clusters = clusters.size();
+    fwrite(&n_clusters, sizeof(n_clusters), 1, fp);
+    fwrite(clusters.data(), clusters.item_size(), clusters.size(), fp);
+}
+
+template <typename ClusterType, typename Enable>
+ClusterVector<ClusterType>
+ClusterFile<ClusterType, Enable>::read_clusters(size_t n_clusters) {
+    if (m_mode != "r") {
+        throw std::runtime_error("File not opened for reading");
+    }
+    if (m_noise_map || m_roi) {
+        return read_clusters_with_cut(n_clusters);
+    } else {
+        return read_clusters_without_cut(n_clusters);
+    }
+}
+
+template <typename ClusterType, typename Enable>
+ClusterVector<ClusterType>
+ClusterFile<ClusterType, Enable>::read_clusters_without_cut(size_t n_clusters) {
+    if (m_mode != "r") {
+        throw std::runtime_error("File not opened for reading");
+    }
+
+    ClusterVector<ClusterType> clusters(n_clusters);
+    clusters.resize(n_clusters);
+
+    int32_t iframe = 0; // frame number needs to be 4 bytes!
+    size_t nph_read = 0;
+    uint32_t nn = m_num_left;
+    uint32_t nph = m_num_left; // number of clusters in frame needs to be 4
+
+    auto buf = clusters.data();
+    // if there are photons left from previous frame read them first
+    if (nph) {
+        if (nph > n_clusters) {
+            // if we have more photons left in the frame then photons to read we
+            // read directly the requested number
+            nn = n_clusters;
+        } else {
+            nn = nph;
+        }
+        nph_read += fread((buf + nph_read), clusters.item_size(), nn, fp);
+        m_num_left = nph - nn; // write back the number of photons left
+    }
+
+    if (nph_read < n_clusters) {
+        // keep on reading frames and photons until reaching n_clusters
+        while (fread(&iframe, sizeof(iframe), 1, fp)) {
+            clusters.set_frame_number(iframe);
+            // read number of clusters in frame
+            if (fread(&nph, sizeof(nph), 1, fp)) {
+                if (nph > (n_clusters - nph_read))
+                    nn = n_clusters - nph_read;
+                else
+                    nn = nph;
+
+                nph_read +=
+                    fread((buf + nph_read), clusters.item_size(), nn, fp);
+                m_num_left = nph - nn;
+            }
+            if (nph_read >= n_clusters)
+                break;
+        }
+    }
+
+    // Resize the vector to the number of clusters.
+    // No new allocation, only change bounds.
+    clusters.resize(nph_read);
+    if (m_gain_map)
+        m_gain_map->apply_gain_map(clusters);
+    return clusters;
+}
+
+template <typename ClusterType, typename Enable>
+ClusterVector<ClusterType>
+ClusterFile<ClusterType, Enable>::read_clusters_with_cut(size_t n_clusters) {
+    ClusterVector<ClusterType> clusters;
+    clusters.resize(n_clusters);
+
+    // if there are photons left from previous frame read them first
+    if (m_num_left) {
+        while (m_num_left && clusters.size() < n_clusters) {
+            ClusterType c = read_one_cluster();
+            if (is_selected(c)) {
+                clusters.push_back(c);
+            }
+        }
+    }
+
+    // we did not have enough clusters left in the previous frame
+    // keep on reading frames until reaching n_clusters
+    if (clusters.size() < n_clusters) {
+        // sanity check
+        if (m_num_left) {
+            throw std::runtime_error(
+                LOCATION + "Entered second loop with clusters left\n");
+        }
+
+        int32_t frame_number = 0; // frame number needs to be 4 bytes!
+        while (fread(&frame_number, sizeof(frame_number), 1, fp)) {
+            if (fread(&m_num_left, sizeof(m_num_left), 1, fp)) {
+                clusters.set_frame_number(
+                    frame_number); // cluster vector will hold the last frame
+                                   // number
+                while (m_num_left && clusters.size() < n_clusters) {
+                    ClusterType c = read_one_cluster();
+                    if (is_selected(c)) {
+                        clusters.push_back(c);
+                    }
+                }
+            }
+
+            // we have enough clusters, break out of the outer while loop
+            if (clusters.size() >= n_clusters)
+                break;
+        }
+    }
+    if (m_gain_map)
+        m_gain_map->apply_gain_map(clusters);
+
+    return clusters;
+}
+
+template <typename ClusterType, typename Enable>
+ClusterType ClusterFile<ClusterType, Enable>::read_one_cluster() {
+    ClusterType c;
+    auto rc = fread(&c, sizeof(c), 1, fp);
+    if (rc != 1) {
+        throw std::runtime_error(LOCATION + "Could not read cluster");
+    }
+    --m_num_left;
+    return c;
+}
+
+template <typename ClusterType, typename Enable>
+ClusterVector<ClusterType> ClusterFile<ClusterType, Enable>::read_frame() {
+    if (m_mode != "r") {
+        throw std::runtime_error(LOCATION + "File not opened for reading");
+    }
+    if (m_noise_map || m_roi) {
+        return read_frame_with_cut();
+    } else {
+        return read_frame_without_cut();
+    }
+}
+
+template <typename ClusterType, typename Enable>
+ClusterVector<ClusterType>
+ClusterFile<ClusterType, Enable>::read_frame_without_cut() {
+    if (m_mode != "r") {
+        throw std::runtime_error("File not opened for reading");
+    }
+    if (m_num_left) {
+        throw std::runtime_error(
+            "There are still photons left in the last frame");
+    }
+    int32_t frame_number;
+    if (fread(&frame_number, sizeof(frame_number), 1, fp) != 1) {
+        throw std::runtime_error(LOCATION + "Could not read frame number");
+    }
+
+    int32_t n_clusters; // Saved as 32bit integer in the cluster file
+    if (fread(&n_clusters, sizeof(n_clusters), 1, fp) != 1) {
+        throw std::runtime_error(LOCATION +
+                                 "Could not read number of clusters");
+    }
+
+    ClusterVector<ClusterType> clusters(n_clusters);
+    clusters.set_frame_number(frame_number);
+
+    if (fread(clusters.data(), clusters.item_size(), n_clusters, fp) !=
+        static_cast<size_t>(n_clusters)) {
+        throw std::runtime_error(LOCATION + "Could not read clusters");
+    }
+    clusters.resize(n_clusters);
+    if (m_gain_map)
+        m_gain_map->apply_gain_map(clusters);
+    return clusters;
+}
+
+template <typename ClusterType, typename Enable>
+ClusterVector<ClusterType>
+ClusterFile<ClusterType, Enable>::read_frame_with_cut() {
+    if (m_mode != "r") {
+        throw std::runtime_error("File not opened for reading");
+    }
+    if (m_num_left) {
+        throw std::runtime_error(
+            "There are still photons left in the last frame");
+    }
+    int32_t frame_number;
+    if (fread(&frame_number, sizeof(frame_number), 1, fp) != 1) {
+        throw std::runtime_error("Could not read frame number");
+    }
+
+    if (fread(&m_num_left, sizeof(m_num_left), 1, fp) != 1) {
+        throw std::runtime_error("Could not read number of clusters");
+    }
+
+    ClusterVector<ClusterType> clusters;
+    clusters.reserve(m_num_left);
+    clusters.set_frame_number(frame_number);
+    while (m_num_left) {
+        ClusterType c = read_one_cluster();
+        if (is_selected(c)) {
+            clusters.push_back(c);
+        }
+    }
+    if (m_gain_map)
+        m_gain_map->apply_gain_map(clusters);
+    return clusters;
+}
+
+template <typename ClusterType, typename Enable>
+bool ClusterFile<ClusterType, Enable>::is_selected(ClusterType &cl) {
+    // Should fail fast
+    if (m_roi) {
+        if (!(m_roi->contains(cl.x, cl.y))) {
+            return false;
+        }
+    }
+
+    auto cluster_size_x = extract_template_arguments<
+        std::remove_reference_t<decltype(cl)>>::cluster_size_x;
+    auto cluster_size_y = extract_template_arguments<
+        std::remove_reference_t<decltype(cl)>>::cluster_size_y;
+
+    size_t cluster_center_index =
+        (cluster_size_x / 2) + (cluster_size_y / 2) * cluster_size_x;
+
+    if (m_noise_map) {
+        auto sum_1x1 = cl.data[cluster_center_index]; // central pixel
+        auto sum_2x2 = cl.max_sum_2x2().first; // highest sum of 2x2 subclusters
+        auto total_sum = cl.sum();             // sum of all pixels
+
+        auto noise =
+            (*m_noise_map)(cl.y, cl.x); // TODO! check if this is correct
+        if (sum_1x1 <= noise || sum_2x2 <= 2 * noise ||
+            total_sum <= 3 * noise) {
+            return false;
+        }
+    }
+    // we passed all checks
+    return true;
+}
+
+} // namespace aare

include/aare/ClusterFileSink.hpp

@@ -0,0 +1,62 @@
+#pragma once
+#include <atomic>
+#include <filesystem>
+#include <thread>
+
+#include "aare/ClusterFinderMT.hpp"
+#include "aare/ClusterVector.hpp"
+#include "aare/ProducerConsumerQueue.hpp"
+
+namespace aare {
+
+template <typename ClusterType,
+          typename = std::enable_if_t<is_cluster_v<ClusterType>>>
+class ClusterFileSink {
+    ProducerConsumerQueue<ClusterVector<ClusterType>> *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<ClusterType> *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<ClusterType, uint16_t, double> *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

@@ -1,15 +1,16 @@
 #pragma once
 #include "aare/core/defs.hpp"
 #include <filesystem>
-#include <string>
 #include <fmt/format.h>
+#include <string>
 
 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);
+        return "frame_number: " + std::to_string(frame_number) +
+               ", n_clusters: " + std::to_string(n_clusters);
     }
 };
 
@@ -24,7 +25,8 @@ struct ClusterV2_ {
                 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) +
+                   ", data: " + data_str;
         }
         return "x: " + std::to_string(x) + ", y: " + std::to_string(y);
     }
@@ -34,27 +36,31 @@ struct ClusterV2 {
     ClusterV2_ cluster;
     int32_t frame_number;
     std::string to_string() const {
-        return "frame_number: " + std::to_string(frame_number) + ", " + cluster.to_string();
+        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
+ * important not: fp always points to the clusters header and does not point to
+ * individual clusters
  *
  */
 class ClusterFileV2 {
-    std::filesystem::path m_fpath; 
+    std::filesystem::path m_fpath;
     std::string m_mode;
     FILE *fp{nullptr};
 
-    void check_open(){
+    void check_open() {
         if (!fp)
-            throw std::runtime_error(fmt::format("File: {} not open", m_fpath.string()));
+            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) {
+    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))
@@ -77,7 +83,7 @@ class ClusterFileV2 {
         check_open();
 
         ClusterHeader header;
-        fread(&header, sizeof(ClusterHeader), 1, fp); 
+        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;
@@ -117,7 +123,7 @@ class ClusterFileV2 {
 
     size_t write(std::vector<std::vector<ClusterV2>> const &clusters) {
         check_open();
-        if (m_mode != "w") 
+        if (m_mode != "w")
             throw std::runtime_error("File not opened in write mode");
 
         size_t n_clusters = 0;

include/aare/ClusterFinder.hpp

@@ -1,216 +1,164 @@
 #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 {
 
-/** enum to define the event types */
-enum eventType {
-    PEDESTAL,            /** pedestal */
-    NEIGHBOUR,           /** neighbour i.e. below threshold, but in the cluster of a photon */
-    PHOTON,              /** photon i.e. above threshold */
-    PHOTON_MAX,          /** maximum of a cluster satisfying the photon conditions */
-    NEGATIVE_PEDESTAL,   /** negative value, will not be accounted for as pedestal in order to
-                             avoid drift of the pedestal towards negative values */
-    UNDEFINED_EVENT = -1 /** undefined */
-};
-
+template <typename ClusterType = Cluster<int32_t, 3, 3>,
+          typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double>
 class ClusterFinder {
+    Shape<2> m_image_size;
+    const PEDESTAL_TYPE m_nSigma;
+    const PEDESTAL_TYPE c2;
+    const PEDESTAL_TYPE c3;
+    Pedestal<PEDESTAL_TYPE> m_pedestal;
+    ClusterVector<ClusterType> m_clusters;
+
+    static const uint8_t ClusterSizeX =
+        extract_template_arguments<ClusterType>::cluster_size_x;
+    static const uint8_t ClusterSizeY =
+        extract_template_arguments<ClusterType>::cluster_size_x;
+    using CT = typename extract_template_arguments<ClusterType>::value_type;
+
   public:
-    ClusterFinder(int cluster_sizeX, int cluster_sizeY, double nSigma = 5.0, double threshold = 0.0)
-        : m_cluster_sizeX(cluster_sizeX), m_cluster_sizeY(cluster_sizeY), m_threshold(threshold), m_nSigma(nSigma) {
+    /**
+     * @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, PEDESTAL_TYPE nSigma = 5.0,
+                  size_t capacity = 1000000)
+        : m_image_size(image_size), m_nSigma(nSigma),
+          c2(sqrt((ClusterSizeY + 1) / 2 * (ClusterSizeX + 1) / 2)),
+          c3(sqrt(ClusterSizeX * ClusterSizeY)),
+          m_pedestal(image_size[0], image_size[1]), m_clusters(capacity) {};
 
-        c2 = sqrt((cluster_sizeY + 1) / 2 * (cluster_sizeX + 1) / 2);
-        c3 = sqrt(cluster_sizeX * cluster_sizeY);
-    };
+    void push_pedestal_frame(NDView<FRAME_TYPE, 2> frame) {
+        m_pedestal.push(frame);
+    }
 
-    template <typename FRAME_TYPE, typename PEDESTAL_TYPE>
-    std::vector<Cluster> find_clusters_without_threshold(NDView<FRAME_TYPE, 2> frame, Pedestal<PEDESTAL_TYPE> &pedestal,
-                                                         bool late_update = false) {
-        struct pedestal_update {
-            int x;
-            int y;
-            FRAME_TYPE value;
-        };
-        std::vector<pedestal_update> pedestal_updates;
+    NDArray<PEDESTAL_TYPE, 2> pedestal() { return m_pedestal.mean(); }
+    NDArray<PEDESTAL_TYPE, 2> noise() { return m_pedestal.std(); }
+    void clear_pedestal() { m_pedestal.clear(); }
 
-        std::vector<Cluster> clusters;
-        std::vector<std::vector<eventType>> eventMask;
-        for (int i = 0; i < frame.shape(0); i++) {
-            eventMask.push_back(std::vector<eventType>(frame.shape(1)));
-        }
-        long double val;
-        long double max;
+    /**
+     * @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<ClusterType>
+    steal_clusters(bool realloc_same_capacity = false) {
+        ClusterVector<ClusterType> tmp = std::move(m_clusters);
+        if (realloc_same_capacity)
+            m_clusters = ClusterVector<ClusterType>(tmp.capacity());
+        else
+            m_clusters = ClusterVector<ClusterType>{};
+        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 = ClusterSizeY / 2;
+        int dx = ClusterSizeX / 2;
+        int has_center_pixel_x =
+            ClusterSizeX %
+            2; // for even sized clusters there is no proper cluster center and
+               // even amount of pixels around the center
+        int has_center_pixel_y = ClusterSizeY % 2;
 
+        m_clusters.set_frame_number(frame_number);
+        std::vector<CT> cluster_data(ClusterSizeX * ClusterSizeY);
         for (int iy = 0; iy < frame.shape(0); iy++) {
             for (int ix = 0; ix < frame.shape(1); ix++) {
-                // initialize max and total
-                max = std::numeric_limits<FRAME_TYPE>::min();
-                long double total = 0;
-                eventMask[iy][ix] = PEDESTAL;
 
-                for (short ir = -(m_cluster_sizeY / 2); ir < (m_cluster_sizeY / 2) + 1; ir++) {
-                    for (short ic = -(m_cluster_sizeX / 2); ic < (m_cluster_sizeX / 2) + 1; ic++) {
-                        if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) {
-                            val = frame(iy + ir, ix + ic) - pedestal.mean(iy + ir, ix + ic);
+                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 + has_center_pixel_y; ir++) {
+                    for (int ic = -dx; ic < dx + has_center_pixel_x; 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;
-                            if (val > max) {
-                                max = val;
-                            }
+                            max = std::max(max, val);
                         }
                     }
                 }
-                auto rms = pedestal.standard_deviation(iy, ix);
-
-                if (frame(iy, ix) - pedestal.mean(iy, ix) < -m_nSigma * rms) {
-                    eventMask[iy][ix] = NEGATIVE_PEDESTAL;
-                    continue;
-                } else if (max > m_nSigma * rms) {
-                    eventMask[iy][ix] = PHOTON;
 
+                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) {
-                    eventMask[iy][ix] = PHOTON;
-                } else{
-                    if (late_update) {
-                        pedestal_updates.push_back({ix, iy, frame(iy, ix)});
-                    } else {
-                        pedestal.push(iy, ix, frame(iy, ix));
-                    }
-                    continue;
+                    // 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
                 }
-                if (eventMask[iy][ix] == PHOTON && (frame(iy, ix) - pedestal.mean(iy, ix)) >= max) {
-                    eventMask[iy][ix] = PHOTON_MAX;
-                    Cluster cluster(m_cluster_sizeX, m_cluster_sizeY, Dtype(typeid(PEDESTAL_TYPE)));
-                    cluster.x = ix;
-                    cluster.y = iy;
-                    short i = 0;
 
-                    for (short ir = -(m_cluster_sizeY / 2); ir < (m_cluster_sizeY / 2) + 1; ir++) {
-                        for (short ic = -(m_cluster_sizeX / 2); ic < (m_cluster_sizeX / 2) + 1; ic++) {
-                            if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) {
-                                PEDESTAL_TYPE tmp = static_cast<PEDESTAL_TYPE>(frame(iy + ir, ix + ic)) -
-                                                    pedestal.mean(iy + ir, ix + ic);
-                                cluster.set<PEDESTAL_TYPE>(i, tmp);
+                // 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 + has_center_pixel_y; ir++) {
+                        for (int ic = -dx; ic < dx + has_center_pixel_y; 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)) -
+                                    static_cast<CT>(
+                                        m_pedestal.mean(iy + ir, ix + ic));
+                                cluster_data[i] =
+                                    tmp; // Watch for out of bounds access
                                 i++;
                             }
                         }
                     }
-                    clusters.push_back(cluster);
+
+                    ClusterType new_cluster{};
+                    new_cluster.x = ix;
+                    new_cluster.y = iy;
+                    std::copy(cluster_data.begin(), cluster_data.end(),
+                              new_cluster.data);
+                    // Add the cluster to the output ClusterVector
+                    m_clusters.push_back(new_cluster);
                 }
             }
         }
-        if (late_update) {
-            for (auto &update : pedestal_updates) {
-                pedestal.push(update.y, update.x, update.value);
-            }
-        }
-        return clusters;
     }
-
-    template <typename FRAME_TYPE, typename PEDESTAL_TYPE>
-    std::vector<Cluster> find_clusters_with_threshold(NDView<FRAME_TYPE, 2> frame, Pedestal<PEDESTAL_TYPE> &pedestal) {
-        assert(m_threshold > 0);
-        std::vector<Cluster> clusters;
-        std::vector<std::vector<eventType>> eventMask;
-        for (int i = 0; i < frame.shape(0); i++) {
-            eventMask.push_back(std::vector<eventType>(frame.shape(1)));
-        }
-        double tthr, tthr1, tthr2;
-
-        NDArray<FRAME_TYPE, 2> rest({frame.shape(0), frame.shape(1)});
-        NDArray<int, 2> nph({frame.shape(0), frame.shape(1)});
-        // convert to n photons
-        // nph = (frame-pedestal.mean()+0.5*m_threshold)/m_threshold; // can be optimized with expression templates?
-        for (int iy = 0; iy < frame.shape(0); iy++) {
-            for (int ix = 0; ix < frame.shape(1); ix++) {
-                auto val = frame(iy, ix) - pedestal.mean(iy, ix);
-                nph(iy, ix) = (val + 0.5 * m_threshold) / m_threshold;
-                nph(iy, ix) = nph(iy, ix) < 0 ? 0 : nph(iy, ix);
-                rest(iy, ix) = val - nph(iy, ix) * m_threshold;
-            }
-        }
-        // iterate over frame pixels
-        for (int iy = 0; iy < frame.shape(0); iy++) {
-            for (int ix = 0; ix < frame.shape(1); ix++) {
-                eventMask[iy][ix] = PEDESTAL;
-                // initialize max and total
-                FRAME_TYPE max = std::numeric_limits<FRAME_TYPE>::min();
-                long double total = 0;
-                if (rest(iy, ix) <= 0.25 * m_threshold) {
-                    pedestal.push(iy, ix, frame(iy, ix));
-                    continue;
-                }
-                eventMask[iy][ix] = NEIGHBOUR;
-                // iterate over cluster pixels around the current pixel (ix,iy)
-                for (short ir = -(m_cluster_sizeY / 2); ir < (m_cluster_sizeY / 2) + 1; ir++) {
-                    for (short ic = -(m_cluster_sizeX / 2); ic < (m_cluster_sizeX / 2) + 1; ic++) {
-                        if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) {
-                            auto val = frame(iy + ir, ix + ic) - pedestal.mean(iy + ir, ix + ic);
-                            total += val;
-                            if (val > max) {
-                                max = val;
-                            }
-                        }
-                    }
-                }
-
-                auto rms = pedestal.standard_deviation(iy, ix);
-                if (m_nSigma == 0) {
-                    tthr = m_threshold;
-                    tthr1 = m_threshold;
-                    tthr2 = m_threshold;
-                } else {
-                    tthr = m_nSigma * rms;
-                    tthr1 = m_nSigma * rms * c3;
-                    tthr2 = m_nSigma * rms * c2;
-
-                    if (m_threshold > 2 * tthr)
-                        tthr = m_threshold - tthr;
-                    if (m_threshold > 2 * tthr1)
-                        tthr1 = tthr - tthr1;
-                    if (m_threshold > 2 * tthr2)
-                        tthr2 = tthr - tthr2;
-                }
-                if (total > tthr1 || max > tthr) {
-                    eventMask[iy][ix] = PHOTON;
-                    nph(iy, ix) += 1;
-                    rest(iy, ix) -= m_threshold;
-                } else {
-                    pedestal.push(iy, ix, frame(iy, ix));
-                    continue;
-                }
-                if (eventMask[iy][ix] == PHOTON && frame(iy, ix) - pedestal.mean(iy, ix) >= max) {
-                    eventMask[iy][ix] = PHOTON_MAX;
-                    Cluster cluster(m_cluster_sizeX, m_cluster_sizeY, Dtype(typeid(FRAME_TYPE)));
-                    cluster.x = ix;
-                    cluster.y = iy;
-                    short i = 0;
-                    for (short ir = -(m_cluster_sizeY / 2); ir < (m_cluster_sizeY / 2) + 1; ir++) {
-                        for (short ic = -(m_cluster_sizeX / 2); ic < (m_cluster_sizeX / 2) + 1; ic++) {
-                            if (ix + ic >= 0 && ix + ic < frame.shape(1) && iy + ir >= 0 && iy + ir < frame.shape(0)) {
-                                auto tmp = frame(iy + ir, ix + ic) - pedestal.mean(iy + ir, ix + ic);
-                                cluster.set<FRAME_TYPE>(i, tmp);
-                                i++;
-                            }
-                        }
-                    }
-                    clusters.push_back(cluster);
-                }
-            }
-        }
-        return clusters;
-    }
-
-  protected:
-    int m_cluster_sizeX;
-    int m_cluster_sizeY;
-    double m_threshold;
-    double m_nSigma;
-    double c2;
-    double c3;
 };
 
 } // namespace aare

include/aare/ClusterFinderMT.hpp

@@ -0,0 +1,274 @@
+#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 ClusterType = Cluster<int32_t, 3, 3>,
+          typename FRAME_TYPE = uint16_t, typename PEDESTAL_TYPE = double>
+class ClusterFinderMT {
+
+    using CT = typename extract_template_arguments<ClusterType>::value_type;
+    size_t m_current_thread{0};
+    size_t m_n_threads{0};
+    using Finder = ClusterFinder<ClusterType, FRAME_TYPE, PEDESTAL_TYPE>;
+    using InputQueue = ProducerConsumerQueue<FrameWrapper>;
+    using OutputQueue = ProducerConsumerQueue<ClusterVector<ClusterType>>;
+    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, 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<ClusterType, FRAME_TYPE, PEDESTAL_TYPE>>(
+                    image_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<ClusterType>> *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,389 @@
+#pragma once
+#include "aare/Cluster.hpp" //TODO maybe store in seperate file !!!
+#include <algorithm>
+#include <array>
+#include <cstddef>
+#include <cstdint>
+#include <numeric>
+#include <vector>
+
+#include <fmt/core.h>
+
+#include "aare/Cluster.hpp"
+#include "aare/NDView.hpp"
+
+namespace aare {
+
+template <typename ClusterType,
+          typename = std::enable_if_t<is_cluster_v<ClusterType>>>
+class ClusterVector; // Forward declaration
+
+/**
+ * @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)
+ */
+#if 0
+template <typename T, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType>
+class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
+
+    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:
+    using value_type = T;
+    using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
+
+    /**
+     * @brief Construct a new ClusterVector object
+     * @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 capacity = 1024, uint64_t frame_number = 0)
+        : m_capacity(capacity), m_frame_number(frame_number) {
+        allocate_buffer(m_capacity);
+    }
+
+    ~ClusterVector() { delete[] m_data; }
+
+    // Move constructor
+    ClusterVector(ClusterVector &&other) noexcept
+        : 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_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
+     */
+    void push_back(const ClusterType &cluster) {
+        if (m_size == m_capacity) {
+            allocate_buffer(m_capacity * 2);
+        }
+        std::byte *ptr = element_ptr(m_size);
+        *reinterpret_cast<CoordType *>(ptr) = cluster.x;
+        ptr += sizeof(CoordType);
+        *reinterpret_cast<CoordType *>(ptr) = cluster.y;
+        ptr += sizeof(CoordType);
+
+        std::memcpy(ptr, cluster.data, ClusterSizeX * ClusterSizeY * sizeof(T));
+
+        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 = ClusterSizeX * ClusterSizeY;
+        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 Sum the pixels in the 2x2 subcluster with the biggest pixel sum in
+     * each cluster
+     * @return std::vector<T> vector of sums for each cluster
+     */ //TODO if underlying container is a vector use std::for_each
+    /*
+    std::vector<T> sum_2x2() {
+        std::vector<T> sums_2x2(m_size);
+
+        for (size_t i = 0; i < m_size; i++) {
+            sums_2x2[i] = at(i).max_sum_2x2;
+        }
+        return sums_2x2;
+    }
+    */
+
+    /**
+     * @brief Return the number of clusters in the vector
+     */
+    size_t size() const { return m_size; }
+
+    uint8_t cluster_size_x() const { return ClusterSizeX; }
+
+    uint8_t cluster_size_y() const { return ClusterSizeY; }
+
+    /**
+     * @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) + ClusterSizeX * ClusterSizeY * 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);
+    }
+
+    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
+     */
+    ClusterType &at(size_t i) {
+        return *reinterpret_cast<ClusterType *>(element_ptr(i));
+    }
+
+    const ClusterType &at(size_t i) const {
+        return *reinterpret_cast<const ClusterType *>(element_ptr(i));
+    }
+
+    template <typename V> const V &at(size_t i) const {
+        return *reinterpret_cast<const 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;
+    }
+};
+#endif
+
+/**
+ * @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, uint8_t ClusterSizeX, uint8_t ClusterSizeY,
+          typename CoordType>
+class ClusterVector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> {
+
+    std::vector<Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>> m_data{};
+    uint64_t m_frame_number{0}; // TODO! Check frame number size and type
+
+  public:
+    using value_type = T;
+    using ClusterType = Cluster<T, ClusterSizeX, ClusterSizeY, CoordType>;
+
+    /**
+     * @brief Construct a new ClusterVector object
+     * @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 capacity = 300, int32_t frame_number = 0)
+        : m_frame_number(frame_number) {
+        m_data.reserve(capacity);
+    }
+
+    // Move constructor
+    ClusterVector(ClusterVector &&other) noexcept
+        : m_data(other.m_data), m_frame_number(other.m_frame_number) {
+        other.m_data.clear();
+    }
+
+    // Move assignment operator
+    ClusterVector &operator=(ClusterVector &&other) noexcept {
+        if (this != &other) {
+            m_data = other.m_data;
+            m_frame_number = other.m_frame_number;
+            other.m_data.clear();
+            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) { m_data.reserve(capacity); }
+
+    void resize(size_t size) { m_data.resize(size); }
+
+    void push_back(const ClusterType &cluster) { m_data.push_back(cluster); }
+
+    ClusterVector &operator+=(const ClusterVector &other) {
+        m_data.insert(m_data.end(), other.begin(), other.end());
+
+        return *this;
+    }
+
+    /**
+     * @brief Return the number of clusters in the vector
+     */
+    size_t size() const { return m_data.size(); }
+
+    uint8_t cluster_size_x() const { return ClusterSizeX; }
+
+    uint8_t cluster_size_y() const { return ClusterSizeY; }
+
+    /**
+     * @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_data.capacity(); }
+
+    const auto begin() const { return m_data.begin(); }
+
+    const auto end() const { return m_data.end(); }
+
+    /**
+     * @brief Return the size in bytes of a single cluster
+     */
+    size_t item_size() const {
+        return 2 * sizeof(CoordType) + ClusterSizeX * ClusterSizeY * sizeof(T);
+    }
+
+    ClusterType *data() { return m_data.data(); }
+    ClusterType const *data() const { return m_data.data(); }
+
+    /**
+     * @brief Return a reference to the i-th cluster casted to type V
+     * @tparam V type of the cluster
+     */
+    ClusterType &at(size_t i) { return m_data[i]; }
+
+    const ClusterType &at(size_t i) const { return m_data[i]; }
+
+    /**
+     * @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;
+    }
+};
+
+} // 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/File.hpp

@@ -1,18 +1,19 @@
 #pragma once
 #include "aare/FileInterface.hpp"
-
+#include <memory>
 
 namespace aare {
 
 /**
- * @brief RAII File class for reading and writing image files in various formats
- * wrapper on a FileInterface to abstract the underlying file format
- * @note documentation for each function is in the FileInterface class
+ * @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 {
-  private:
-    FileInterface *file_impl;
-    bool is_npy;
+    std::unique_ptr<FileInterface> file_impl;
 
   public:
     /**
@@ -24,37 +25,42 @@ class File {
      * @throws std::invalid_argument if the file mode is not supported
      *
      */
-    File(const std::filesystem::path &fname, const std::string &mode, const FileConfig &cfg = {});
-    void write(Frame &frame, sls_detector_header header = {});
-    Frame read();
-    Frame iread(size_t frame_number);
-    std::vector<Frame> read(size_t n_frames);
+    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(size_t frame_index);
-    size_t bytes_per_frame();
-    size_t pixels_per_frame();
-    void seek(size_t frame_number);
-    size_t tell() const;
+    
+    size_t 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;
-    size_t bitdepth() const;
-    size_t bytes_per_pixel() const;
-    void set_total_frames(size_t total_frames);
+
     DetectorType detector_type() const;
-    xy geometry() const;
 
-    /**
-     * @brief Move constructor
-     * @param other File object to move from
-     */
-    File(File &&other) noexcept;
-
-    /**
-     * @brief destructor: will only delete the FileInterface object
-     */
-    ~File();
+    
 };
 
 } // namespace aare

include/aare/FileInterface.hpp

@@ -33,7 +33,7 @@ struct FileConfig {
     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) +
+               ", 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) + " }";
     }
@@ -47,32 +47,24 @@ struct FileConfig {
 class FileInterface {
   public:
     /**
-     * @brief write a frame to the file
-     * @param frame frame to write
-     * @return void
-     * @throws std::runtime_error if the function is not implemented
-     */
-    // virtual void write(Frame &frame) = 0;
-
-    /**
-     * @brief write a vector of frames to the file
-     * @param frames vector of frames to write
-     * @return void
-     */
-    // virtual void write(std::vector<Frame> &frames) = 0;
-
-    /**
-     * @brief read one frame from the file at the current position
+     * @brief  one frame from the file at the current position
      * @return Frame
      */
-    virtual Frame read() = 0;
+    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(size_t n_frames) = 0; // Is this the right interface?
+    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
@@ -142,55 +134,28 @@ class FileInterface {
      */
     virtual size_t bitdepth() const = 0;
 
-    /**
-     * @brief read one frame from the file at the given frame number
-     * @param frame_number frame number to read
-     * @return frame
-     */
-    Frame iread(size_t frame_number) {
-        auto old_pos = tell();
-        seek(frame_number);
-        Frame tmp = read();
-        seek(old_pos);
-        return tmp;
-    };
 
-    /**
-     * @brief read n_frames from the file starting at the given frame number
-     * @param frame_number frame number to start reading from
-     * @param n_frames number of frames to read
-     * @return vector of frames
-     */
-    std::vector<Frame> iread(size_t frame_number, size_t n_frames) {
-        auto old_pos = tell();
-        seek(frame_number);
-        std::vector<Frame> tmp = read(n_frames);
-        seek(old_pos);
-        return tmp;
-    }
-    DetectorType detector_type() const { return m_type; }
+    virtual DetectorType detector_type() const = 0;
 
     // function to query the data type of the file
     /*virtual DataType dtype = 0; */
 
     virtual ~FileInterface() = default;
 
-    void set_total_frames(size_t total_frames) { m_total_frames = total_frames; }
-
   protected:
     std::string m_mode{};
-    std::filesystem::path m_fname{};
-    std::filesystem::path m_base_path{};
-    std::string m_base_name{}, m_ext{};
-    int m_findex{};
-    size_t m_total_frames{};
-    size_t max_frames_per_file{};
-    std::string version{};
-    DetectorType m_type{DetectorType::Unknown};
-    size_t m_rows{};
-    size_t m_cols{};
-    size_t m_bitdepth{};
-    size_t current_frame{};
+    // 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/FilePtr.hpp

@@ -0,0 +1,30 @@
+#pragma once
+#include <cstdio>
+#include <filesystem>
+
+namespace aare {
+
+/**
+ * \brief RAII wrapper for FILE pointer
+ */
+class FilePtr {
+    FILE *fp_{nullptr};
+
+  public:
+    FilePtr() = default;
+    FilePtr(const std::filesystem::path& fname, const std::string& mode);
+    FilePtr(const FilePtr &) = delete;            // we don't want a copy
+    FilePtr &operator=(const FilePtr &) = delete; // since we handle a resource
+    FilePtr(FilePtr &&other);
+    FilePtr &operator=(FilePtr &&other);
+    FILE *get();
+    int64_t tell();
+    void seek(int64_t offset, int whence = SEEK_SET) {
+        if (fseek(fp_, offset, whence) != 0)
+            throw std::runtime_error("Error seeking in file");
+    }
+    std::string error_msg();
+    ~FilePtr();
+};
+
+} // 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

@@ -11,31 +11,49 @@
 namespace aare {
 
 /**
- * @brief Frame class to represent a single frame of data
- * model class
- * should be able to work with streams coming from files or network
+ * @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);
-    Frame(const std::byte *bytes, uint32_t rows, uint32_t cols, Dtype dtype);
-    ~Frame() noexcept;
 
-    // disable copy and assignment
-    Frame &operator=(const Frame &other)=delete;
-    Frame(const Frame &other)=delete;
+    /**
+     * @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;
 
-    // explicit copy
-    Frame copy() const;
+
+    Frame clone() const; //<- Explicit copy
 
     uint32_t rows() const;
     uint32_t cols() const;
@@ -45,32 +63,62 @@ class Frame {
     size_t bytes() const;
     std::byte *data() const;
 
-    std::byte *get(uint32_t row, uint32_t col);
+    /**
+     * @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;
 
-    // TODO! can we, or even want to remove the template?
+    /**
+     * @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());
+        std::memcpy(m_data + (row * m_cols + col) * m_dtype.bytes(), &data,
+                    m_dtype.bytes());
     }
-    template <typename T> T get_t(uint32_t row, uint32_t col) {
+    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());
+        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)};
+        std::array<int64_t, 2> shape = {static_cast<int64_t>(m_rows),
+                                        static_cast<int64_t>(m_cols)};
         T *data = reinterpret_cast<T *>(m_data);
         return NDView<T, 2>(data, shape);
     }
 
-    template <typename T> NDArray<T> image() { return NDArray<T>(this->view<T>()); }
+    /** 
+     * @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/GainMap.hpp

@@ -0,0 +1,58 @@
+/************************************************
+ * @file ApplyGainMap.hpp
+ * @short function to apply gain map of image size to a vector of clusters
+ ***********************************************/
+
+#pragma once
+#include "aare/Cluster.hpp"
+#include "aare/ClusterVector.hpp"
+#include "aare/NDArray.hpp"
+#include "aare/NDView.hpp"
+#include <memory>
+
+namespace aare {
+
+class GainMap {
+
+  public:
+    explicit GainMap(const NDArray<double, 2> &gain_map)
+        : m_gain_map(gain_map) {};
+
+    explicit GainMap(const NDView<double, 2> gain_map) {
+        m_gain_map = NDArray<double, 2>(gain_map);
+    }
+
+    template <typename ClusterType,
+              typename = std::enable_if_t<is_cluster_v<ClusterType>>>
+    void apply_gain_map(ClusterVector<ClusterType> &clustervec) {
+        // in principle we need to know the size of the image for this lookup
+        size_t ClusterSizeX = clustervec.cluster_size_x();
+        size_t ClusterSizeY = clustervec.cluster_size_y();
+
+        using T = typename ClusterVector<ClusterType>::value_type;
+
+        int64_t index_cluster_center_x = ClusterSizeX / 2;
+        int64_t index_cluster_center_y = ClusterSizeY / 2;
+        for (size_t i = 0; i < clustervec.size(); i++) {
+            auto &cl = clustervec.at(i);
+
+            if (cl.x > 0 && cl.y > 0 && cl.x < m_gain_map.shape(1) - 1 &&
+                cl.y < m_gain_map.shape(0) - 1) {
+                for (size_t j = 0; j < ClusterSizeX * ClusterSizeY; j++) {
+                    size_t x = cl.x + j % ClusterSizeX - index_cluster_center_x;
+                    size_t y = cl.y + j / ClusterSizeX - index_cluster_center_y;
+                    cl.data[j] = cl.data[j] * static_cast<T>(m_gain_map(y, x));
+                }
+            } else {
+                memset(cl.data, 0,
+                       ClusterSizeX * ClusterSizeY *
+                           sizeof(T)); // clear edge clusters
+            }
+        }
+    }
+
+  private:
+    NDArray<double, 2> m_gain_map{};
+};
+
+} // end of namespace aare

include/aare/Interpolator.hpp

@@ -0,0 +1,132 @@
+#pragma once
+
+#include "aare/CalculateEta.hpp"
+#include "aare/Cluster.hpp"
+#include "aare/ClusterFile.hpp" //Cluster_3x3
+#include "aare/ClusterVector.hpp"
+#include "aare/NDArray.hpp"
+#include "aare/NDView.hpp"
+#include "aare/algorithm.hpp"
+
+namespace aare {
+
+struct Photon {
+    double x;
+    double y;
+    double energy;
+};
+
+class Interpolator {
+    NDArray<double, 3> m_ietax;
+    NDArray<double, 3> m_ietay;
+
+    NDArray<double, 1> m_etabinsx;
+    NDArray<double, 1> m_etabinsy;
+    NDArray<double, 1> m_energy_bins;
+
+  public:
+    Interpolator(NDView<double, 3> etacube, NDView<double, 1> xbins,
+                 NDView<double, 1> ybins, NDView<double, 1> ebins);
+    NDArray<double, 3> get_ietax() { return m_ietax; }
+    NDArray<double, 3> get_ietay() { return m_ietay; }
+
+    template <typename ClusterType,
+              typename Eanble = std::enable_if_t<is_cluster_v<ClusterType>>>
+    std::vector<Photon> interpolate(const ClusterVector<ClusterType> &clusters);
+};
+
+// TODO: generalize to support any clustertype!!! otherwise add std::enable_if_t
+// to only take Cluster2x2 and Cluster3x3
+template <typename ClusterType, typename Enable>
+std::vector<Photon>
+Interpolator::interpolate(const ClusterVector<ClusterType> &clusters) {
+    std::vector<Photon> photons;
+    photons.reserve(clusters.size());
+
+    if (clusters.cluster_size_x() == 3 || clusters.cluster_size_y() == 3) {
+        for (size_t i = 0; i < clusters.size(); i++) {
+
+            auto cluster = clusters.at(i);
+            auto eta = calculate_eta2(cluster);
+
+            Photon photon;
+            photon.x = cluster.x;
+            photon.y = cluster.y;
+            photon.energy = eta.sum;
+
+            // auto ie = nearest_index(m_energy_bins, photon.energy)-1;
+            // auto ix = nearest_index(m_etabinsx, eta.x)-1;
+            // auto iy = nearest_index(m_etabinsy, eta.y)-1;
+            // Finding the index of the last element that is smaller
+            // should work fine as long as we have many bins
+            auto ie = last_smaller(m_energy_bins, photon.energy);
+            auto ix = last_smaller(m_etabinsx, eta.x);
+            auto iy = last_smaller(m_etabinsy, eta.y);
+
+            // fmt::print("ex: {}, ix: {}, iy: {}\n", ie, ix, iy);
+
+            double dX, dY;
+            // cBottomLeft = 0,
+            // cBottomRight = 1,
+            // cTopLeft = 2,
+            // cTopRight = 3
+            switch (eta.c) {
+            case cTopLeft:
+                dX = -1.;
+                dY = 0;
+                break;
+            case cTopRight:;
+                dX = 0;
+                dY = 0;
+                break;
+            case cBottomLeft:
+                dX = -1.;
+                dY = -1.;
+                break;
+            case cBottomRight:
+                dX = 0.;
+                dY = -1.;
+                break;
+            }
+            photon.x += m_ietax(ix, iy, ie) * 2 + dX;
+            photon.y += m_ietay(ix, iy, ie) * 2 + dY;
+            photons.push_back(photon);
+        }
+    } else if (clusters.cluster_size_x() == 2 ||
+               clusters.cluster_size_y() == 2) {
+        for (size_t i = 0; i < clusters.size(); i++) {
+            auto cluster = clusters.at(i);
+            auto eta = calculate_eta2(cluster);
+
+            Photon photon;
+            photon.x = cluster.x;
+            photon.y = cluster.y;
+            photon.energy = eta.sum;
+
+            // Now do some actual interpolation.
+            // Find which energy bin the cluster is in
+            //  auto ie = nearest_index(m_energy_bins, photon.energy)-1;
+            //  auto ix = nearest_index(m_etabinsx, eta.x)-1;
+            //  auto iy = nearest_index(m_etabinsy, eta.y)-1;
+            // Finding the index of the last element that is smaller
+            // should work fine as long as we have many bins
+            auto ie = last_smaller(m_energy_bins, photon.energy);
+            auto ix = last_smaller(m_etabinsx, eta.x);
+            auto iy = last_smaller(m_etabinsy, eta.y);
+
+            photon.x += m_ietax(ix, iy, ie) *
+                        2; // eta goes between 0 and 1 but we could move the hit
+                           // anywhere in the 2x2
+            photon.y += m_ietay(ix, iy, ie) * 2;
+            photons.push_back(photon);
+        }
+
+    } else {
+        throw std::runtime_error(
+            "Only 3x3 and 2x2 clusters are supported for interpolation");
+    }
+
+    return photons;
+}
+
+} // namespace aare

include/aare/JungfrauDataFile.hpp

@@ -0,0 +1,106 @@
+#pragma once
+#include <cstdint>
+#include <filesystem>
+#include <vector>
+
+#include "aare/FilePtr.hpp"
+#include "aare/defs.hpp"
+#include "aare/NDArray.hpp"
+#include "aare/FileInterface.hpp"
+namespace aare {
+
+
+struct JungfrauDataHeader{
+    uint64_t framenum;
+    uint64_t bunchid;
+};
+
+class JungfrauDataFile : public FileInterface {
+
+    size_t m_rows{};     //!< number of rows in the image, from find_frame_size();
+    size_t m_cols{};     //!< number of columns in the image, from find_frame_size();
+    size_t m_bytes_per_frame{}; //!< number of bytes per frame excluding header
+    size_t m_total_frames{};    //!< total number of frames in the series of files
+    size_t m_offset{}; //!< file index of the first file, allow starting at non zero file
+    size_t m_current_file_index{};   //!< The index of the open file
+    size_t m_current_frame_index{};  //!< The index of the current frame (with reference to all files)
+
+    std::vector<size_t> m_last_frame_in_file{}; //!< Used for seeking to the correct file
+    std::filesystem::path m_path; //!< path to the files
+    std::string m_base_name;      //!< base name used for formatting file names
+
+    FilePtr m_fp; //!< RAII wrapper for a FILE* 
+
+    
+    using pixel_type = uint16_t;
+    static constexpr size_t header_size = sizeof(JungfrauDataHeader); 
+    static constexpr size_t n_digits_in_file_index = 6; //!< to format file names
+
+  public:
+    JungfrauDataFile(const std::filesystem::path &fname);
+
+    std::string base_name() const;      //!< get the base name of the file (without path and extension)
+    size_t bytes_per_frame() override;           
+    size_t pixels_per_frame() override;    
+    size_t bytes_per_pixel() const;  
+    size_t bitdepth() const override;
+    void seek(size_t frame_index) override;  //!< seek to the given frame index (note not byte offset)
+    size_t tell() override;            //!< get the frame index of the file pointer
+    size_t total_frames() const override;
+    size_t rows() const override;
+    size_t cols() const override;
+    std::array<ssize_t,2> shape() const;
+    size_t n_files() const; //!< get the number of files in the series. 
+
+    // Extra functions needed for FileInterface
+    Frame read_frame() override;
+    Frame read_frame(size_t frame_number) override;
+    std::vector<Frame> read_n(size_t n_frames=0) override;
+    void read_into(std::byte *image_buf) override;
+    void read_into(std::byte *image_buf, size_t n_frames) override;
+    size_t frame_number(size_t frame_index) override;
+    DetectorType detector_type() const override;
+
+    /**
+     * @brief Read a single frame from the file into the given buffer.  
+     * @param image_buf buffer to read the frame into. (Note the caller is responsible for allocating the buffer)
+     * @param header pointer to a JungfrauDataHeader or nullptr to skip header)
+     */
+    void read_into(std::byte *image_buf, JungfrauDataHeader *header = nullptr);
+
+    /**
+     * @brief Read a multiple frames from the file into the given buffer.  
+     * @param image_buf buffer to read the frame into. (Note the caller is responsible for allocating the buffer)
+     * @param n_frames number of frames to read
+     * @param header pointer to a JungfrauDataHeader or nullptr to skip header)
+     */
+    void read_into(std::byte *image_buf, size_t n_frames, JungfrauDataHeader *header = nullptr);
+    
+    /** 
+     * @brief Read a single frame from the file into the given NDArray
+     * @param image NDArray to read the frame into.
+     */
+    void read_into(NDArray<uint16_t>* image, JungfrauDataHeader* header = nullptr);
+
+    JungfrauDataHeader read_header();
+    std::filesystem::path current_file() const { return fpath(m_current_file_index+m_offset); }
+
+    
+    private:
+    /**
+     * @brief Find the size of the frame in the file. (256x256, 256x1024, 512x1024)
+     * @param fname path to the file
+     * @throws std::runtime_error if the file is empty or the size cannot be determined
+     */
+     void find_frame_size(const std::filesystem::path &fname);
+
+
+     void parse_fname(const std::filesystem::path &fname);
+     void scan_files();
+     void open_file(size_t file_index);
+     std::filesystem::path fpath(size_t frame_index) const;
+
+
+  };
+
+} // namespace aare

include/aare/NDArray.hpp

@@ -7,6 +7,7 @@ memory.
 TODO! Add expression templates for operators
 
 */
+#include "aare/ArrayExpr.hpp"
 #include "aare/NDView.hpp"
 
 #include <algorithm>
@@ -20,36 +21,80 @@ TODO! Add expression templates for operators
 
 namespace aare {
 
-template <typename T, int64_t Ndim = 2> class NDArray {
-  public:
-    NDArray() : shape_(), strides_(c_strides<Ndim>(shape_)), data_(nullptr){};
 
+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_]){};
+          size_(std::accumulate(shape_.begin(), shape_.end(), 1,
+                                std::multiplies<>())),
+          data_(new T[size_]) {}
 
-    NDArray(std::array<int64_t, Ndim> shape, T value) : NDArray(shape) { this->operator=(value); }
 
-    /* When constructing from a NDView we need to copy the data since
-    NDArray expect to own its data, and span is just a view*/
-    explicit NDArray(NDView<T, Ndim> span) : NDArray(span.shape()) {
-        std::copy(span.begin(), span.end(), begin());
-        // fmt::print("NDArray(NDView<T, Ndim> span)\n");
+    /**
+     * @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_) {
+        : shape_(other.shape_), strides_(c_strides<Ndim>(shape_)),
+          size_(other.size_), data_(other.data_) {
+        other.reset(); // TODO! is this necessary?
 
-        other.reset();
-        // fmt::print("NDArray(NDArray &&other)\n");
     }
 
     // Copy constructor
     NDArray(const NDArray &other)
-        : shape_(other.shape_), strides_(c_strides<Ndim>(shape_)), size_(other.size_), data_(new T[size_]) {
+        : shape_(other.shape_), strides_(c_strides<Ndim>(shape_)),
+          size_(other.size_), data_(new T[size_]) {
         std::copy(other.data_, other.data_ + size_, data_);
-        // fmt::print("NDArray(const NDArray &other)\n");
+    }
+
+    // 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_; }
@@ -57,19 +102,33 @@ template <typename T, int64_t Ndim = 2> class NDArray {
     auto begin() { return data_; }
     auto end() { return data_ + size_; }
 
+    auto begin() const { return data_; }
+    auto end() const { return data_ + size_; }
+
     using value_type = T;
 
     NDArray &operator=(NDArray &&other) noexcept; // Move assign
     NDArray &operator=(const NDArray &other);     // Copy assign
-
-    NDArray operator+(const NDArray &other);
     NDArray &operator+=(const NDArray &other);
-    NDArray operator-(const NDArray &other);
     NDArray &operator-=(const NDArray &other);
-    NDArray operator*(const NDArray &other);
     NDArray &operator*=(const NDArray &other);
-    NDArray operator/(const NDArray &other);
+
+    //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()) {
@@ -98,6 +157,11 @@ template <typename T, int64_t Ndim = 2> class NDArray {
 
     NDArray &operator&=(const T & /*mask*/);
 
+  
+
+
+    
+
     void sqrt() {
         for (int i = 0; i < size_; ++i) {
             data_[i] = std::sqrt(data_[i]);
@@ -106,38 +170,50 @@ template <typename T, int64_t Ndim = 2> class NDArray {
 
     NDArray &operator++(); // pre inc
 
-    template <typename... Ix> std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) {
+    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 {
+    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) {
+    template <typename... Ix>
+    std::enable_if_t<sizeof...(Ix) == Ndim, T> value(Ix... index) {
         return data_[element_offset(strides_, index...)];
     }
 
-    T &operator()(int i) { return data_[i]; }
-    const T &operator()(int i) const { return data_[i]; }
+    // 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_); }
-    uint64_t size() const { return size_; }
+    ssize_t size() const { return static_cast<ssize_t>(size_); }
     size_t total_bytes() const { return size_ * sizeof(T); }
     std::array<int64_t, Ndim> shape() const noexcept { return shape_; }
     int64_t shape(int64_t i) const noexcept { return shape_[i]; }
     std::array<int64_t, Ndim> strides() const noexcept { return strides_; }
     size_t bitdepth() const noexcept { return sizeof(T) * 8; }
+
     std::array<int64_t, Ndim> byte_strides() const noexcept {
         auto byte_strides = strides_;
         for (auto &val : byte_strides)
             val *= sizeof(T);
         return byte_strides;
-        // return strides_;
     }
 
-    NDView<T, Ndim> span() const { return NDView<T, Ndim>{data_, shape_}; }
+    /**
+     * @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();
@@ -149,16 +225,12 @@ template <typename T, int64_t Ndim = 2> class NDArray {
         std::fill(shape_.begin(), shape_.end(), 0);
         std::fill(strides_.begin(), strides_.end(), 0);
     }
-
-  private:
-    std::array<int64_t, Ndim> shape_;
-    std::array<int64_t, Ndim> strides_;
-    uint64_t size_{};
-    T *data_;
 };
 
 // Move assign
-template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator=(NDArray<T, Ndim> &&other) noexcept {
+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_;
@@ -170,15 +242,11 @@ template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator
     return *this;
 }
 
-template <typename T, int64_t Ndim> NDArray<T, Ndim> NDArray<T, Ndim>::operator+(const NDArray &other) {
-    NDArray result(*this);
-    result += other;
-    return result;
-}
-template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const NDArray<T, Ndim> &other) {
+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) {
+        for (size_t i = 0; i < size_; ++i) {
             data_[i] += other.data_[i];
         }
         return *this;
@@ -186,13 +254,8 @@ template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator
     throw(std::runtime_error("Shape of ImageDatas must match"));
 }
 
-template <typename T, int64_t Ndim> NDArray<T, Ndim> NDArray<T, Ndim>::operator-(const NDArray &other) {
-    NDArray result{*this};
-    result -= other;
-    return result;
-}
-
-template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator-=(const NDArray<T, Ndim> &other) {
+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) {
@@ -202,13 +265,9 @@ template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator
     }
     throw(std::runtime_error("Shape of ImageDatas must match"));
 }
-template <typename T, int64_t Ndim> NDArray<T, Ndim> NDArray<T, Ndim>::operator*(const NDArray &other) {
-    NDArray result = *this;
-    result *= other;
-    return result;
-}
 
-template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator*=(const NDArray<T, Ndim> &other) {
+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) {
@@ -219,36 +278,17 @@ template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator
     throw(std::runtime_error("Shape of ImageDatas must match"));
 }
 
-template <typename T, int64_t Ndim> NDArray<T, Ndim> NDArray<T, Ndim>::operator/(const NDArray &other) {
-    NDArray result = *this;
-    result /= other;
-    return result;
-}
-
-template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator&=(const T &mask) {
+template <typename T, int64_t Ndim>
+NDArray<T, Ndim> &NDArray<T, Ndim>::operator&=(const T &mask) {
     for (auto it = begin(); it != end(); ++it)
         *it &= mask;
     return *this;
 }
 
-// template <typename T, int64_t Ndim>
-// NDArray<T, Ndim>& NDArray<T, Ndim>::operator/=(const NDArray<T, Ndim>&
-// other)
-// {
-//     //check shape
-//     if (shape_ == other.shape_) {
-//         for (int i = 0; i < size_; ++i) {
-//             data_[i] /= other.data_[i];
-//         }
-//         return *this;
-//     } else {
-//         throw(std::runtime_error("Shape of ImageDatas must match"));
-//     }
-// }
-
-template <typename T, int64_t Ndim> NDArray<bool, Ndim> NDArray<T, Ndim>::operator>(const NDArray &other) {
+template <typename T, int64_t Ndim>
+NDArray<bool, Ndim> NDArray<T, Ndim>::operator>(const NDArray &other) {
     if (shape_ == other.shape_) {
-        NDArray<bool> result{shape_};
+        NDArray<bool, Ndim> result{shape_};
         for (int i = 0; i < size_; ++i) {
             result(i) = (data_[i] > other.data_[i]);
         }
@@ -257,7 +297,8 @@ template <typename T, int64_t Ndim> NDArray<bool, Ndim> NDArray<T, Ndim>::operat
     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) {
+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_;
@@ -269,7 +310,8 @@ template <typename T, int64_t Ndim> NDArray<T, Ndim> &NDArray<T, Ndim>::operator
     return *this;
 }
 
-template <typename T, int64_t Ndim> bool NDArray<T, Ndim>::operator==(const NDArray<T, Ndim> &other) const {
+template <typename T, int64_t Ndim>
+bool NDArray<T, Ndim>::operator==(const NDArray<T, Ndim> &other) const {
     if (shape_ != other.shape_)
         return false;
 
@@ -280,67 +322,94 @@ template <typename T, int64_t Ndim> bool NDArray<T, Ndim>::operator==(const NDAr
     return true;
 }
 
-template <typename T, int64_t Ndim> bool NDArray<T, Ndim>::operator!=(const NDArray<T, Ndim> &other) const {
+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++() {
+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) {
+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) {
+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) {
+
+
+
+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) {
+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) {
+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) {
+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) {
+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) {
+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) {
+template <typename T, int64_t Ndim>
+NDArray<T, Ndim> NDArray<T, Ndim>::operator*(const T &value) {
     NDArray result = *this;
     result *= value;
     return result;
 }
-template <typename T, int64_t Ndim> void NDArray<T, Ndim>::Print() {
-    if (shape_[0] < 20 && shape_[1] < 20)
-        Print_all();
-    else
-        Print_some();
+// template <typename T, int64_t Ndim> void NDArray<T, Ndim>::Print() {
+//     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) {
@@ -360,7 +429,8 @@ template <typename T, int64_t Ndim> void NDArray<T, Ndim>::Print_some() {
     }
 }
 
-template <typename T, int64_t Ndim> void save(NDArray<T, Ndim> &img, std::string &pathname) {
+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));
@@ -368,7 +438,8 @@ template <typename T, int64_t Ndim> void save(NDArray<T, Ndim> &img, std::string
 }
 
 template <typename T, int64_t Ndim>
-NDArray<T, Ndim> load(const std::string &pathname, std::array<int64_t, Ndim> shape) {
+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);
@@ -377,4 +448,6 @@ NDArray<T, Ndim> load(const std::string &pathname, std::array<int64_t, Ndim> sha
     return img;
 }
 
+
+
 } // namespace aare

include/aare/NDView.hpp

@@ -1,4 +1,7 @@
 #pragma once
+#include "aare/defs.hpp"
+#include "aare/ArrayExpr.hpp"
+
 #include <algorithm>
 #include <array>
 #include <cassert>
@@ -45,7 +48,7 @@ template <int64_t Ndim> std::array<int64_t, Ndim> make_array(const std::vector<i
     return arr;
 }
 
-template <typename T, int64_t Ndim = 2> class NDView {
+template <typename T, int64_t Ndim = 2> class NDView : public ArrayExpr<NDView<T, Ndim>, Ndim> {
   public:
     NDView() = default;
     ~NDView() = default;
@@ -68,12 +71,14 @@ template <typename T, int64_t Ndim = 2> class NDView {
         return buffer_[element_offset(strides_, index...)];
     }
 
-    uint64_t size() const { return size_; }
+    ssize_t size() const { return static_cast<ssize_t>(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]; }
 
@@ -94,6 +99,15 @@ template <typename T, int64_t Ndim = 2> class NDView {
 
     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() != static_cast<ssize_t>(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;
@@ -121,7 +135,7 @@ template <typename T, int64_t Ndim = 2> class NDView {
         return *this;
     }
 
-    auto &shape() { return shape_; }
+    auto &shape() const { return shape_; }
     auto shape(int64_t i) const { return shape_[i]; }
 
     T *data() { return buffer_; }
@@ -156,4 +170,18 @@ template <typename T, int64_t Ndim> void NDView<T, Ndim>::print_all() const {
     }
 }
 
+
+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

@@ -31,9 +31,10 @@ class NumpyFile : public FileInterface {
     explicit NumpyFile(const std::filesystem::path &fname, const std::string &mode = "r", FileConfig cfg = {});
 
     void write(Frame &frame);
-    Frame read() override { return get_frame(this->current_frame++); }
+    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(size_t n_frames) override;
+    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; };
@@ -46,6 +47,8 @@ class NumpyFile : public FileInterface {
     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
@@ -103,6 +106,10 @@ class NumpyFile : public FileInterface {
     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);

include/aare/Pedestal.hpp

@@ -6,116 +6,204 @@
 
 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_freeze(false), m_samples(n_samples),           m_cur_samples(NDArray<uint32_t, 2>({rows, cols}, 0)),m_sum(NDArray<SUM_TYPE, 2>({rows, cols})),
- m_sum2(NDArray<SUM_TYPE, 2>({rows, cols})) {
+        : 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() {
-        NDArray<SUM_TYPE, 2> mean_array({m_rows, m_cols});
-        for (uint32_t i = 0; i < m_rows * m_cols; i++) {
-            mean_array(i / m_cols, i % m_cols) = mean(i / m_cols, i % m_cols);
+        return 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 mean_array;
+        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);
+            variance_array(i / m_cols, i % m_cols) =
+                variance(i / m_cols, i % m_cols);
         }
         return variance_array;
     }
 
-    NDArray<SUM_TYPE, 2> standard_deviation() {
+    
+
+    NDArray<SUM_TYPE, 2> 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) = standard_deviation(i / m_cols, i % m_cols);
+            standard_deviation_array(i / m_cols, i % m_cols) =
+                std(i / m_cols, i % m_cols);
         }
 
         return standard_deviation_array;
     }
+
+    
+
     void clear() {
-        for (uint32_t i = 0; i < m_rows * m_cols; i++) {
-            clear(i / m_cols, i % m_cols);
-        }
+        m_sum = 0;
+        m_sum2 = 0;
+        m_cur_samples = 0;
+        m_mean = 0;
     }
 
-    /*
-     * index level operations
-     */
-    SUM_TYPE mean(const uint32_t row, const uint32_t col) const {
-        if (m_cur_samples(row, col) == 0) {
-            return 0.0;
-        }
-        return m_sum(row, col) / m_cur_samples(row, col);
-    }
-    SUM_TYPE variance(const uint32_t row, const uint32_t col) const {
-        if (m_cur_samples(row, col) == 0) {
-            return 0.0;
-        }
-        return m_sum2(row, col) / m_cur_samples(row, col) - mean(row, col) * mean(row, col);
-    }
-    SUM_TYPE standard_deviation(const uint32_t row, const uint32_t col) const { return std::sqrt(variance(row, col)); }
+    
 
     void clear(const uint32_t row, const uint32_t col) {
         m_sum(row, col) = 0;
         m_sum2(row, col) = 0;
         m_cur_samples(row, col) = 0;
+        m_mean(row, col) = 0;
     }
-    // frame level operations
+    
+
+
     template <typename T> void push(NDView<T, 2> frame) {
         assert(frame.size() == m_rows * m_cols);
-        // TODO: test the effect of #pragma omp parallel for
-        for (uint32_t index = 0; index < m_rows * m_cols; index++) {
-            push<T>(index / m_cols, index % m_cols, frame(index));
+
+        // 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));
+        assert(frame.rows() == static_cast<size_t>(m_rows) &&
+               frame.cols() == static_cast<size_t>(m_cols));
         push<T>(frame.view<T>());
     }
 
     // getter functions
-    inline uint32_t rows() const { return m_rows; }
-    inline uint32_t cols() const { return m_cols; }
-    inline uint32_t n_samples() const { return m_samples; }
-    inline NDArray<uint32_t, 2> cur_samples() const { return m_cur_samples; }
-    inline NDArray<SUM_TYPE, 2> get_sum() const { return m_sum; }
-    inline NDArray<SUM_TYPE, 2> get_sum2() const { return m_sum2; }
+    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> inline void push(const uint32_t row, const uint32_t col, const T val_) {
-        if (m_freeze) {
-            return;
-        }
+    // 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_);
-        const uint32_t idx = index(row, col);
-        if (m_cur_samples(idx) < m_samples) {
-            m_sum(idx) += val;
-            m_sum2(idx) += val * val;
-            m_cur_samples(idx)++;
+        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(idx) += val - m_sum(idx) / m_cur_samples(idx);
-            m_sum2(idx) += val * val - m_sum2(idx) / m_cur_samples(idx);
+            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);
         }
     }
-    inline uint32_t index(const uint32_t row, const uint32_t col) const { return row * m_cols + col; };
-    void set_freeze(bool freeze) { m_freeze = freeze; }
 
-  private:
-    uint32_t m_rows;
-    uint32_t m_cols;
-    bool m_freeze;
-    uint32_t m_samples;
-    NDArray<uint32_t, 2> m_cur_samples;
-    NDArray<SUM_TYPE, 2> m_sum;
-    NDArray<SUM_TYPE, 2> m_sum2;
+    /**
+     * @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

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

include/aare/SubFile.hpp

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

include/aare/VarClusterFinder.hpp

@@ -7,7 +7,7 @@
 
 #include "aare/NDArray.hpp"
 
-const int MAX_CLUSTER_SIZE = 200;
+const int MAX_CLUSTER_SIZE = 50;
 namespace aare {
 
 template <typename T> class VarClusterFinder {
@@ -49,7 +49,7 @@ template <typename T> class VarClusterFinder {
     int check_neighbours(int i, int j);
 
   public:
-    VarClusterFinder(image_shape shape, T threshold)
+    VarClusterFinder(Shape<2> shape, T threshold)
         : shape_(shape), labeled_(shape, 0), peripheral_labeled_(shape, 0), binary_(shape), threshold_(threshold) {
         hits.reserve(2000);
     }
@@ -226,7 +226,7 @@ template <typename T> void VarClusterFinder<T>::single_pass(NDView<T, 2> img) {
 
 template <typename T> void VarClusterFinder<T>::first_pass() {
 
-    for (int i = 0; i < original_.size(); ++i) {
+    for (ssize_t i = 0; i < original_.size(); ++i) {
         if (use_noise_map)
             threshold_ = 5 * noiseMap(i);
         binary_(i) = (original_(i) > threshold_);
@@ -250,17 +250,17 @@ template <typename T> void VarClusterFinder<T>::first_pass() {
 
 template <typename T> void VarClusterFinder<T>::second_pass() {
 
-    for (int64_t i = 0; i != labeled_.size(); ++i) {
-        auto current_label = labeled_(i);
-        if (current_label != 0) {
-            auto it = child.find(current_label);
+    for (ssize_t i = 0; i != labeled_.size(); ++i) {
+        auto cl = labeled_(i);
+        if (cl != 0) {
+            auto it = child.find(cl);
             while (it != child.end()) {
-                current_label = it->second;
-                it = child.find(current_label);
+                cl = it->second;
+                it = child.find(cl);
                 // do this once before doing the second pass?
                 // all values point to the final one...
             }
-            labeled_(i) = current_label;
+            labeled_(i) = cl;
         }
     }
 }
@@ -271,7 +271,7 @@ template <typename T> void VarClusterFinder<T>::store_clusters() {
     // Do we always have monotonic increasing
     // labels? Then vector?
     // here the translation is label -> Hit
-    std::unordered_map<int, Hit> h_size;
+    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
@@ -280,7 +280,7 @@ template <typename T> void VarClusterFinder<T>::store_clusters() {
                 // (i+1 < shape_[0] and labeled_(i+1, j) != 0) or
                 // (j+1 < shape_[1] and labeled_(i, j+1) != 0)
             ) {
-                Hit &record = h_size[labeled_(i, j)];
+                Hit &record = h_map[labeled_(i, j)];
                 if (record.size < MAX_CLUSTER_SIZE) {
                     record.rows[record.size] = i;
                     record.cols[record.size] = j;
@@ -300,7 +300,7 @@ template <typename T> void VarClusterFinder<T>::store_clusters() {
         }
     }
 
-    for (const auto &h : h_size)
+    for (const auto &h : h_map)
         hits.push_back(h.second);
 }
 

include/aare/algorithm.hpp

@@ -0,0 +1,111 @@
+
+#pragma once
+#include <algorithm>
+#include <array>
+#include <vector>
+#include <aare/NDArray.hpp>
+
+namespace aare {
+/**
+ * @brief Index of the last element that is smaller than val.
+ * Requires a sorted array. Uses >= for ordering. If all elements
+ * are smaller it returns the last element and if all elements are
+ * larger it returns the first element.
+ * @param first iterator to the first element
+ * @param last iterator to the last element
+ * @param val value to compare
+ * @return index of the last element that is smaller than val
+ * 
+ */
+template <typename T>
+size_t last_smaller(const T* first, const T* last, T val) {
+    for (auto iter = first+1; iter != last; ++iter) {
+        if (*iter >= val) {
+            return std::distance(first, iter-1);
+        }
+    }
+    return std::distance(first, last-1);
+}
+
+template <typename T>
+size_t last_smaller(const NDArray<T, 1>& arr, T val) {
+    return last_smaller(arr.begin(), arr.end(), val);
+}
+
+template <typename T>
+size_t last_smaller(const std::vector<T>& vec, T val) {
+    return last_smaller(vec.data(), vec.data()+vec.size(), val);
+}
+
+/**
+ * @brief Index of the first element that is larger than val.
+ * Requires a sorted array. Uses > for ordering. If all elements
+ * are larger it returns the first element and if all elements are
+ * smaller it returns the last element.
+ * @param first iterator to the first element
+ * @param last iterator to the last element
+ * @param val value to compare
+ * @return index of the first element that is larger than val
+ */
+template <typename T>
+size_t first_larger(const T* first, const T* last, T val) {
+    for (auto iter = first; iter != last; ++iter) {
+        if (*iter > val) {
+            return std::distance(first, iter);
+        }
+    }
+    return std::distance(first, last-1);
+}
+
+template <typename T>
+size_t first_larger(const NDArray<T, 1>& arr, T val) {
+    return first_larger(arr.begin(), arr.end(), val);
+}
+
+template <typename T>
+size_t first_larger(const std::vector<T>& vec, T val) {
+    return first_larger(vec.data(), vec.data()+vec.size(), val);
+}
+
+/**
+ * @brief Index of the nearest element to val.
+ * Requires a sorted array. If there is no difference it takes the first element.
+ * @param first iterator to the first element
+ * @param last iterator to the last element
+ * @param val value to compare
+ * @return index of the nearest element
+ */
+template <typename T>
+size_t nearest_index(const T* first, const T* last, T val) {
+    auto iter = std::min_element(first, last,
+    [val](T a, T b) {
+        return std::abs(a - val) < std::abs(b - val);
+    });
+    return std::distance(first, iter);
+}
+
+template <typename T>
+size_t nearest_index(const NDArray<T, 1>& arr, T val) {
+    return nearest_index(arr.begin(), arr.end(), val);
+}
+
+template <typename T>
+size_t nearest_index(const std::vector<T>& vec, T val) {
+    return nearest_index(vec.data(), vec.data()+vec.size(), val);
+}
+
+template <typename T, size_t N>
+size_t nearest_index(const std::array<T,N>& arr, T val) {
+    return nearest_index(arr.data(), arr.data()+arr.size(), val);
+}
+
+template <typename T>
+std::vector<T> cumsum(const std::vector<T>& vec) {
+    std::vector<T> result(vec.size());
+    std::partial_sum(vec.begin(), vec.end(), result.begin());
+    return result;
+}
+
+
+
+} // 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

@@ -1,11 +1,9 @@
 #pragma once
 
 #include "aare/Dtype.hpp"
-// #include "aare/utils/logger.hpp"
 
 #include <array>
 #include <stdexcept>
-
 #include <cassert>
 #include <cstdint>
 #include <cstring>
@@ -14,70 +12,101 @@
 #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__) + ":"
+#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 {
 
-class Cluster {
+inline constexpr size_t bits_per_byte = 8;
+
+void assert_failed(const std::string &msg);
+
+
+
+class DynamicCluster {
   public:
     int cluster_sizeX;
     int cluster_sizeY;
     int16_t x;
     int16_t y;
-    Dtype dt;
+    Dtype dt; // 4 bytes
 
   private:
     std::byte *m_data;
 
   public:
-    Cluster(int cluster_sizeX_, int cluster_sizeY_, Dtype dt_ = Dtype(typeid(int32_t)))
-        : cluster_sizeX(cluster_sizeX_), cluster_sizeY(cluster_sizeY_), dt(dt_) {
+    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()]{};
     }
-    Cluster() : Cluster(3, 3) {}
-    Cluster(const Cluster &other) : Cluster(other.cluster_sizeX, other.cluster_sizeY, other.dt) {
+    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());
     }
-    Cluster &operator=(const Cluster &other) {
+    DynamicCluster &operator=(const DynamicCluster &other) {
         if (this == &other)
             return *this;
-        this->~Cluster();
-        new (this) Cluster(other);
+        this->~DynamicCluster();
+        new (this) DynamicCluster(other);
         return *this;
     }
-    Cluster(Cluster &&other) noexcept
-        : cluster_sizeX(other.cluster_sizeX), cluster_sizeY(other.cluster_sizeY), x(other.x), y(other.y), dt(other.dt),
-          m_data(other.m_data) {
+    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);
     }
-    ~Cluster() { delete[] m_data; }
+    ~DynamicCluster() { delete[] m_data; }
     template <typename T> T get(int idx) {
-        (sizeof(T) == dt.bytes()) ? 0 : throw std::invalid_argument("[ERROR] Type size mismatch");
+        (sizeof(T) == dt.bytes())
+            ? 0
+            : throw std::invalid_argument("[ERROR] Type size mismatch");
         return *reinterpret_cast<T *>(m_data + idx * dt.bytes());
     }
     template <typename T> auto set(int idx, T val) {
-        (sizeof(T) == dt.bytes()) ? 0 : throw std::invalid_argument("[ERROR] Type size mismatch");
-        return memcpy(m_data + idx * dt.bytes(), &val, (size_t)dt.bytes());
+        (sizeof(T) == dt.bytes())
+            ? 0
+            : throw std::invalid_argument("[ERROR] Type size mismatch");
+        return memcpy(m_data + idx * dt.bytes(), &val, dt.bytes());
     }
-    // auto x() const { return x; }
-    // auto y() const { return y; }
-    // auto x(int16_t x_) { return x = x_; }
-    // auto y(int16_t y_) { return y = y_; }
 
     template <typename T> std::string to_string() const {
-        (sizeof(T) == dt.bytes()) ? 0 : throw std::invalid_argument("[ERROR] Type size mismatch");
-        std::string s = "x: " + std::to_string(x) + " y: " + std::to_string(y) + "\nm_data: [";
+        (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 += std::to_string(
+                     *reinterpret_cast<T *>(m_data + i * dt.bytes())) +
+                 " ";
         }
         s += "]";
         return s;
@@ -85,17 +114,19 @@ class Cluster {
     /**
      * @brief size of the cluster in bytes when saved to a file
      */
-    size_t size() const { return cluster_sizeX * cluster_sizeY ; }
+    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(); }
+    auto end() const {
+        return m_data + cluster_sizeX * cluster_sizeY * dt.bytes();
+    }
     std::byte *data() { return m_data; }
 };
 
 /**
  * @brief header contained in parts of frames
  */
-struct sls_detector_header {
+struct DetectorHeader {
     uint64_t frameNumber;
     uint32_t expLength;
     uint32_t packetNumber;
@@ -117,40 +148,119 @@ struct sls_detector_header {
         }
         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";
+        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 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) + " }"; }
+    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; }
+    bool contains(int64_t x, int64_t y) const {
+        return x >= xmin && x < xmax && y >= ymin && y < ymax;
+    }
+  };
+
+
 using dynamic_shape = std::vector<int64_t>;
 
-enum class DetectorType { Jungfrau, Eiger, Mythen3, Moench, ChipTestBoard, Unknown };
+//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 <class T> std::string ToString(T arg) { return T(arg); }
 
 template <> DetectorType StringTo(const std::string & /*name*/);
-template <> std::string toString(DetectorType arg);
+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/json.hpp

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

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/libzmq_cmake_version.patch

@@ -0,0 +1,18 @@
+diff --git a/CMakeLists.txt b/CMakeLists.txt
+index dd3d8eb9..c0187747 100644
+--- a/CMakeLists.txt
++++ b/CMakeLists.txt
+@@ -1,11 +1,8 @@
+ # CMake build script for ZeroMQ
+ project(ZeroMQ)
+ 
+-if(${CMAKE_SYSTEM_NAME} STREQUAL Darwin)
+-  cmake_minimum_required(VERSION 3.0.2)
+-else()
+-  cmake_minimum_required(VERSION 2.8.12)
+-endif()
++cmake_minimum_required(VERSION 3.15)
++message(STATUS "Patched cmake version")
+ 
+ include(CheckIncludeFiles)
+ include(CheckCCompilerFlag)

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,23 @@
+[build-system]
+requires = ["scikit-build-core>=0.10", "pybind11", "numpy"]
+build-backend = "scikit_build_core.build"
+
+[project]
+name = "aare"
+version = "2025.4.1"
+
+
+
+
+[tool.scikit-build]
+cmake.verbose = true
+
+[tool.scikit-build.cmake.define]
+AARE_PYTHON_BINDINGS = "ON"
+AARE_SYSTEM_LIBRARIES = "ON"
+AARE_INSTALL_PYTHONEXT = "ON"
+
+[tool.pytest.ini_options]
+markers = [
+    "files: marks tests that need additional data (deselect with '-m \"not files\"')",
+]

python/CMakeLists.txt

@@ -1,34 +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.11.0
+        GIT_TAG        v2.13.0
     )
     FetchContent_MakeAvailable(pybind11)
 else()
-    find_package(pybind11 2.11 REQUIRED)
+    find_package(pybind11 2.13 REQUIRED)
 endif()
 
-
+# Add the compiled python extension
 pybind11_add_module(
-    _aare 
-    src/module.cpp
+    _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}  )
+    configure_file(${FILE} ${CMAKE_BINARY_DIR}/${FILE}  )
 endforeach(FILE ${PYTHON_FILES})
 
+set_target_properties(_aare PROPERTIES
+    LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/aare
+)
 
-configure_file( examples/play.py  
-    ${CMAKE_BINARY_DIR}/play.py 
-)
+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

@@ -1,3 +1,24 @@
+# Make the compiled classes that live in _aare available from aare.
+from . import _aare
 
 
-from _aare import File
+from ._aare import File, RawMasterFile, RawSubFile, JungfrauDataFile
+from ._aare import Pedestal_d, Pedestal_f, ClusterFinder_Cluster3x3i, VarClusterFinder
+from ._aare import DetectorType
+from ._aare import ClusterFile_Cluster3x3i as ClusterFile
+from ._aare import hitmap
+from ._aare import ROI
+
+# from ._aare import ClusterFinderMT, ClusterCollector, ClusterFileSink, ClusterVector_i
+
+from ._aare import fit_gaus, fit_pol1
+from ._aare import Interpolator
+from .CtbRawFile import CtbRawFile
+from .RawFile import RawFile
+from .ScanParameters import ScanParameters
+
+from .utils import random_pixels, random_pixel, flat_list, add_colorbar
+
+
+#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,36 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1 import make_axes_locatable
+
+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]
+
+def add_colorbar(ax, im, size="5%", pad=0.05):
+    """Add a colorbar with the same height as the image."""
+    divider = make_axes_locatable(ax)
+    cax = divider.append_axes("right", size=size, pad=pad)
+    plt.colorbar(im, cax=cax)
+    return ax, im, cax

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

@@ -1,14 +1,79 @@
-import matplotlib.pyplot as plt
+import sys
+sys.path.append('/home/l_msdetect/erik/aare/build')
+
+from aare._aare import ClusterVector_i, Interpolator
+
+import pickle 
 import numpy as np
-plt.ion()
+import matplotlib.pyplot as plt
+import boost_histogram as bh
+import torch
+import math
+import time
 
-import aare
-from pathlib import Path
 
-p = Path('/Users/erik/data/aare_test_data/jungfrau/jungfrau_single_master_0.json')
 
-f = aare.File(p)
-frame = f.read_frame()
+def gaussian_2d(mx, my, sigma = 1, res=100, grid_size = 2):
+    """
+    Generate a 2D gaussian as position mx, my, with sigma=sigma. 
+    The gaussian is placed on a 2x2 pixel matrix with resolution 
+    res in one dimesion.
+    """
+    x = torch.linspace(0, pixel_size*grid_size, res)
+    x,y = torch.meshgrid(x,x, indexing="ij")
+    return 1 / (2*math.pi*sigma**2) * \
+      torch.exp(-((x - my)**2 / (2*sigma**2) + (y - mx)**2 / (2*sigma**2)))
 
-fig, ax = plt.subplots()
-im = ax.imshow(frame, cmap='viridis')
+scale = 1000 #Scale factor when converting to integer
+pixel_size = 25 #um
+grid = 2
+resolution = 100
+sigma_um = 10
+xa = np.linspace(0,grid*pixel_size,resolution)
+ticks = [0, 25, 50]
+
+hit = np.array((20,20))
+etahist_fname = "/home/l_msdetect/erik/tmp/test_hist.pkl"
+
+local_resolution = 99
+grid_size = 3
+xaxis = np.linspace(0,grid_size*pixel_size, local_resolution)
+t = gaussian_2d(hit[0],hit[1], grid_size = grid_size, sigma = 10, res = local_resolution)
+pixels = t.reshape(grid_size, t.shape[0] // grid_size, grid_size, t.shape[1] // grid_size).sum(axis = 3).sum(axis = 1)
+pixels = pixels.numpy()
+pixels = (pixels*scale).astype(np.int32)
+v = ClusterVector_i(3,3)
+v.push_back(1,1, pixels)
+
+with open(etahist_fname, "rb") as f:
+        hist = pickle.load(f)
+eta = hist.view().copy()
+etabinsx = np.array(hist.axes.edges.T[0].flat)
+etabinsy = np.array(hist.axes.edges.T[1].flat)
+ebins = np.array(hist.axes.edges.T[2].flat)
+p = Interpolator(eta, etabinsx[0:-1], etabinsy[0:-1], ebins[0:-1])
+
+
+
+
+#Generate the hit
+
+
+
+
+tmp = p.interpolate(v)
+print(f'tmp:{tmp}')
+pos = np.array((tmp['x'], tmp['y']))*25
+
+
+print(pixels)
+fig, ax = plt.subplots(figsize = (7,7))
+ax.pcolormesh(xaxis, xaxis, t)
+ax.plot(*pos, 'o')
+ax.set_xticks([0,25,50,75])
+ax.set_yticks([0,25,50,75])
+ax.set_xlim(0,75)
+ax.set_ylim(0,75)
+ax.grid()
+print(f'{hit=}')
+print(f'{pos=}')