133 Commits

Author SHA1 Message Date
7d6223d52d Cluster finder improvements (#113) 2024-12-16 14:42:18 +01:00
da67f58323 Cluster finder improvements (#112) 2024-12-16 14:26:35 +01:00
e6098c02ef bumped version 2024-12-16 14:24:46 +01:00
29b1dc8df3 missing header 2024-12-13 14:57:36 +01:00
f88b53387f WIP 2024-12-12 17:58:04 +01:00
a0f481c0ee mod pedestal 2024-12-12 14:34:10 +01:00
b3a9e9576b WIP 2024-12-11 16:27:36 +01:00
60534add92 WIP 2024-12-11 09:54:33 +01:00
7f2a23d5b1 accumulating clusters in one array 2024-12-10 22:00:12 +01:00
6a150e8d98 WIP 2024-12-10 17:21:05 +01:00
b43003966f build pkg on all branched deploy docs on main 2024-11-29 16:41:42 +01:00
c2d039a5bd fix conda build 2024-11-29 16:37:42 +01:00
6fd52f6b8d added missing enums (#111)
- Missing enums
- Matching values to slsDetectorPackage
- tests
2024-11-29 15:28:19 +01:00
659f1f36c5 AARE_INSTALL_PYTHONEXT (#109)
- added AARE_INSTALL_PYTHONEXT option to install also python files in
aare folder
2024-11-29 15:28:02 +01:00
0047d15de1 removed print flip 2024-11-29 15:00:18 +01:00
a1b7fb8fc8 added missing enums 2024-11-29 14:56:39 +01:00
2e4a491d7a CMAKE_INSTALL_PREFIX not needed to specify destination folder and removed 2024-11-29 14:38:32 +01:00
fdce2f69b9 python 3.10 required in cmake 2024-11-29 11:07:05 +01:00
ada4d41f4a bugfix on iteration and returning master file (#110) 2024-11-29 08:52:04 +01:00
115dfc0abf bugfix on iteration and returning master file 2024-11-28 21:14:40 +01:00
31b834c3fd added AARE_INSTALL_PYTHONEXT option to install python in make install, which also installs the python files in the aare folder 2024-11-28 15:18:13 +01:00
feed4860a6 Developer (#108)
- Support for very old moench
- read_n in RawFile
2024-11-27 21:22:01 +01:00
0df8e4bb7d added support for old old moench files 2024-11-27 16:27:55 +01:00
8bf9ac55ce modified read_n also for File and RawFile 2024-11-27 09:31:57 +01:00
2d33fd4813 Streamlined build, new transforms (#106) 2024-11-26 16:27:00 +01:00
996a8861f6 roll back conda-build 2024-11-26 15:53:06 +01:00
06670a7e24 read_n returns remaining frames (#105)
Modified read_n to return the number of frames available if less than
the number of frames requested.

```python
#f is a CtbRawFile containing 10 frames

f.read_n(7) # you get 7 frames
f.read_n(7) # you get 3 frames
f.read_n(7) # RuntimeError
```

Also added support for chunk_size when iterating over a file:

```python
# The file contains 10 frames


with CtbRawFile(fname, chunk_size = 7) as f:
    for headers, frames in f:
        #do something with the data
        # 1 iteration 7 frames
        # 2 iteration 3 frames
        # 3 iteration stops
```
2024-11-26 14:07:21 +01:00
8e3d997bed read_n returns remaining frames 2024-11-26 12:07:17 +01:00
a3f813f9b4 Modified moench05 transform (#103)
Moench05 transforms: 
- moench05: Works with the updated firmware and better data compression
(adcenable10g=0xFF0F)
- moench05_old: Works with the previous data and can be used with
adcenable10g=0xFFFFFFFF
- moench05_1g: For the 1g data acquisition only with adcenable=0x2202
2024-11-26 09:02:33 +01:00
d48482e9da Modified moench05 transform: new firmware (moench05), legacy firmware (moench05_old), 1g readout (moench05_1g) 2024-11-25 16:39:08 +01:00
8f729fc83e Developer (#102) 2024-11-21 10:27:26 +01:00
f9a2d49244 removed extra print 2024-11-21 10:22:22 +01:00
9f7cdbcb48 conversion warnings 2024-11-18 18:18:55 +01:00
3b0e13e41f added links (#101) 2024-11-18 16:19:15 +01:00
3af8182998 added links 2024-11-18 16:18:29 +01:00
99e829fd06 Latest changes (#100) 2024-11-18 15:42:37 +01:00
47e867fc1a Merge branch 'main' into developer 2024-11-18 15:38:15 +01:00
8ea4372cf1 fix 2024-11-18 15:33:38 +01:00
75f83e5e3b detecting need to link with stdfs 2024-11-18 15:33:09 +01:00
30d05f9203 detecting need to link with stdfs 2024-11-18 15:19:57 +01:00
37d3dfcf71 WIP 2024-11-18 14:46:28 +01:00
35c6706b3c docs 2024-11-18 14:39:46 +01:00
9ab61cac4e deps in pkg 2024-11-18 11:47:26 +01:00
13394c3a61 cmake targets 2024-11-18 11:30:33 +01:00
088288787a Merge branch 'developer' of github.com:slsdetectorgroup/aare into developer 2024-11-18 09:22:36 +01:00
9d4459eb8c linking json with PUBLIC to avoid errors 2024-11-18 09:22:28 +01:00
95ff77c8fc Cluster reading, expression templates (#99)
Co-authored-by: froejdh_e <erik.frojdh@psi.ch>
2024-11-15 16:32:36 +01:00
62a14dda13 Merge branch 'main' into developer 2024-11-15 16:19:34 +01:00
632c2ee0c8 bumped version 2024-11-15 16:15:04 +01:00
17f8d28019 frame reading for cluster file 2024-11-15 16:13:46 +01:00
e77b615293 Added expression templates (#98)
- Works with NDArray
- Works with NDView
2024-11-15 15:17:52 +01:00
0d058274d5 WIP 2024-11-14 17:03:16 +01:00
5cde7a99b5 WIP 2024-11-14 17:02:48 +01:00
dcedb4fb13 added missing header 2024-11-14 16:37:24 +01:00
7ffd732d98 ported reading clusters (#95) 2024-11-14 16:22:38 +01:00
fbaf9dce89 Developer (#94) 2024-11-14 08:03:18 +01:00
dc889dab76 removed subfile from cmake 2024-11-14 07:48:59 +01:00
cb94d079af Merge branch 'developer' of github.com:slsdetectorgroup/aare into developer 2024-11-14 07:42:00 +01:00
13b2cb40b6 docs and reorder 2024-11-14 07:41:50 +01:00
17917ac7ea Merge branch 'main' into developer 2024-11-12 16:44:15 +01:00
db936b6357 improved documentation 2024-11-12 16:39:03 +01:00
2ee1a5583e WIP 2024-11-12 09:27:01 +01:00
349e3af8e1 Brining in changes (#93) 2024-11-11 19:59:55 +01:00
a0b6c4cc03 Merge branch 'main' into developer 2024-11-11 18:52:23 +01:00
5f21759c8c removed prints, bumped version 2024-11-11 18:22:18 +01:00
ecf1b2a90b WIP 2024-11-11 17:13:48 +01:00
b172c7aa0a starting work on ROI 2024-11-07 16:24:48 +01:00
d8d1f0c517 Taking v1 as the first release (#92)
- file reading
- decoding master file
2024-11-07 10:14:20 +01:00
d5fb823ae4 added numpy variants 2024-11-07 09:16:49 +01:00
9c220bff51 added simple decoding of scan parameters 2024-11-07 08:14:33 +01:00
b2e5c71f9c MH02 1-4 counters 2024-11-06 21:32:03 +01:00
cbfd1f0b6c ClusterFinder 2024-11-06 12:41:41 +01:00
5b2809d6b0 working Moench03 detector type 2024-11-06 10:13:56 +01:00
4bb8487e2c added moench03 back 2024-11-06 09:18:57 +01:00
1cc7690f9a discard partial 2024-11-06 09:13:40 +01:00
25812cb291 RawFile is now using RawMasterFile 2024-11-06 09:10:09 +01:00
654c1db3f4 WIP 2024-11-05 16:01:22 +01:00
2efb763242 func to prop 2024-11-05 16:00:11 +01:00
7f244e22a2 extra methods in CtbRawFile 2024-11-05 15:55:17 +01:00
d98b45235f optional 2024-11-05 14:37:35 +01:00
80a39415de added CtbRawFile 2024-11-05 14:36:18 +01:00
b8a4498379 WIP 2024-10-31 18:03:17 +01:00
49da039ff9 working on 05 2024-10-31 15:35:43 +01:00
563c39c0dd decoding of old Moench03 2024-10-31 11:53:24 +01:00
ae1166b908 WIP 2024-10-31 10:39:52 +01:00
ec61132296 WIP 2024-10-31 10:29:07 +01:00
cee0d71b9c added check to prevent segfault on missing subfile 2024-10-31 09:43:41 +01:00
19c6a4091f improved docs and added PixelMap 2024-10-31 08:56:12 +01:00
92d9c28c73 numpy in conda env for docs 2024-10-30 18:35:54 +01:00
b7e6962e44 added numpy as dep 2024-10-30 18:09:29 +01:00
13ac6b0f37 added missing numpy dependency 2024-10-30 17:54:34 +01:00
79d924c2a3 docs and version bump 2024-10-30 17:48:50 +01:00
9b733fd0ec WIP 2024-10-30 17:41:45 +01:00
6505f37d87 added type bindings 2024-10-30 17:40:47 +01:00
a466887064 added variable cluster finder 2024-10-30 17:22:57 +01:00
dde92b993f xml back in 2024-10-30 16:29:43 +01:00
1b61155c5c another try 2024-10-30 16:15:44 +01:00
738934f2a0 added github token 2024-10-30 16:12:08 +01:00
6b8f2478b6 added deploy 2024-10-30 16:04:11 +01:00
41fbddb750 pinned sphinx version 2024-10-30 15:56:27 +01:00
504e8b4565 updated doxyfile 2024-10-30 15:53:40 +01:00
acdcaac338 fmt 2024-10-30 15:39:04 +01:00
8b43011fa1 modified action 2024-10-30 15:37:23 +01:00
801adccbd7 updated path for docs 2024-10-30 15:28:07 +01:00
da5ba034b8 WIP 2024-10-30 15:18:48 +01:00
1cbded04f8 doxygen 2024-10-30 15:14:55 +01:00
9b33ad0ee8 pybind 2024-10-30 15:13:20 +01:00
1f539a234b forgot json 2024-10-30 15:11:52 +01:00
29a42507d7 WIP 2024-10-30 15:09:35 +01:00
5035c20aa4 added action for docs 2024-10-30 15:07:34 +01:00
f754e0f769 file reading 2024-10-30 14:53:50 +01:00
be019b9769 updated readme 2024-10-30 10:26:53 +01:00
af4f000fe7 fetch content for json 2024-10-30 09:36:41 +01:00
b37f4845cf cmake defaults 2024-10-30 08:58:42 +01:00
b037aebc5f update 2024-10-30 08:36:38 +01:00
dea5aaf9cf slight mod 2024-10-30 08:33:22 +01:00
4cc6aa9d40 updated workflows 2024-10-30 08:11:38 +01:00
c3a5d22f83 added anaconda-client 2024-10-29 17:54:02 +01:00
a8afa04129 updated workflow 2024-10-29 17:50:44 +01:00
eb855fb9a3 updated workflow 2024-10-29 17:49:21 +01:00
b4fe044679 WIP 2024-10-29 17:00:58 +01:00
082d793161 WIP 2024-10-29 16:46:02 +01:00
9f29f173ff updated path 2024-10-29 13:09:14 +01:00
8a10bcbbdb workflow 2024-10-29 13:07:53 +01:00
c509e29b52 building with scikit build 2024-10-29 11:19:20 +01:00
1a16d4522e WIP 2024-10-28 16:50:38 +01:00
8a435cbe9b WIP 2024-10-28 16:26:14 +01:00
7f9151f270 WIP 2024-10-28 13:37:58 +01:00
abb1d20ca3 WIP 2024-10-28 12:25:47 +01:00
a4fb217e3f Files and structure for python interface 2024-10-28 11:22:12 +01:00
5d643dc133 added cluster finder 2024-10-25 16:18:36 +02:00
54dd88f070 added documentation 2024-10-25 13:54:36 +02:00
b1b020ad60 WIP 2024-10-25 10:23:34 +02:00
110 changed files with 12352 additions and 0 deletions

7
.clang-format Normal file
View File

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

68
.github/workflows/build_docs.yml vendored Normal file
View File

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

44
.github/workflows/deploy.yml vendored Normal file
View File

@ -0,0 +1,44 @@
name: Build pkgs and deploy if on main
on:
push:
branches:
- main
- 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: Enable upload
if: github.ref == 'refs/heads/main'
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}

24
.gitignore vendored Normal file
View File

@ -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__/*

439
CMakeLists.txt Normal file
View File

@ -0,0 +1,439 @@
cmake_minimum_required(VERSION 3.14)
project(aare
VERSION 1.0.0
DESCRIPTION "Data processing library for PSI detectors"
HOMEPAGE_URL "https://github.com/slsdetectorgroup/aare"
LANGUAGES C CXX
)
set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED ON)
set(CMAKE_CXX_EXTENSIONS OFF)
if (${CMAKE_VERSION} VERSION_GREATER "3.24")
cmake_policy(SET CMP0135 NEW) #Fetch content download timestamp
endif()
cmake_policy(SET CMP0079 NEW)
include(GNUInstallDirs)
include(FetchContent)
#Set default build type if none was specified
include(cmake/helpers.cmake)
default_build_type("Release")
set_std_fs_lib()
message(STATUS "Extra linking to fs lib:${STD_FS_LIB}")
set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH})
# General options
option(AARE_PYTHON_BINDINGS "Build python bindings" ON)
option(AARE_TESTS "Build tests" OFF)
option(AARE_BENCHMARKS "Build benchmarks" OFF)
option(AARE_EXAMPLES "Build examples" OFF)
option(AARE_IN_GITHUB_ACTIONS "Running in Github Actions" OFF)
option(AARE_DOCS "Build documentation" OFF)
option(AARE_VERBOSE "Verbose output" OFF)
option(AARE_CUSTOM_ASSERT "Use custom assert" OFF)
option(AARE_INSTALL_PYTHONEXT "Install the python extension in the install tree under CMAKE_INSTALL_PREFIX/aare/" OFF)
option(AARE_ASAN "Enable AddressSanitizer" OFF)
# Configure which of the dependencies to use FetchContent for
option(AARE_FETCH_FMT "Use FetchContent to download fmt" ON)
option(AARE_FETCH_PYBIND11 "Use FetchContent to download pybind11" ON)
option(AARE_FETCH_CATCH "Use FetchContent to download catch2" ON)
option(AARE_FETCH_JSON "Use FetchContent to download nlohmann::json" ON)
option(AARE_FETCH_ZMQ "Use FetchContent to download libzmq" ON)
#Convenience option to use system libraries only (no FetchContent)
option(AARE_SYSTEM_LIBRARIES "Use system libraries" OFF)
if(AARE_SYSTEM_LIBRARIES)
message(STATUS "Build using system libraries")
set(AARE_FETCH_FMT OFF CACHE BOOL "Disabled FetchContent for FMT" FORCE)
set(AARE_FETCH_PYBIND11 OFF CACHE BOOL "Disabled FetchContent for pybind11" FORCE)
set(AARE_FETCH_CATCH OFF CACHE BOOL "Disabled FetchContent for catch2" FORCE)
set(AARE_FETCH_JSON OFF CACHE BOOL "Disabled FetchContent for nlohmann::json" FORCE)
set(AARE_FETCH_ZMQ OFF CACHE BOOL "Disabled FetchContent for libzmq" FORCE)
endif()
if(AARE_VERBOSE)
add_compile_definitions(AARE_VERBOSE)
endif()
if(AARE_CUSTOM_ASSERT)
add_compile_definitions(AARE_CUSTOM_ASSERT)
endif()
if(AARE_BENCHMARKS)
add_subdirectory(benchmarks)
endif()
set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
if(AARE_FETCH_ZMQ)
# Fetchcontent_Declare is deprecated need to find a way to update this
# for now setting the policy to old is enough
if (${CMAKE_VERSION} VERSION_GREATER_EQUAL "3.30")
cmake_policy(SET CMP0169 OLD)
endif()
FetchContent_Declare(
libzmq
GIT_REPOSITORY https://github.com/zeromq/libzmq.git
GIT_TAG v4.3.4
)
# Disable unwanted options from libzmq
set(BUILD_TESTS OFF CACHE BOOL "Switch off libzmq test build")
set(BUILD_SHARED OFF CACHE BOOL "Switch off libzmq shared libs")
set(WITH_PERF_TOOL OFF CACHE BOOL "")
set(ENABLE_CPACK OFF CACHE BOOL "")
set(ENABLE_CLANG OFF CACHE BOOL "")
set(ENABLE_CURVE OFF CACHE BOOL "")
set(ENABLE_DRAFTS OFF CACHE BOOL "")
# TODO! Verify that this is what we want to do in aare
# Using GetProperties and Populate to be able to exclude zmq
# from install (not possible with FetchContent_MakeAvailable(libzmq))
FetchContent_GetProperties(libzmq)
if(NOT libzmq_POPULATED)
FetchContent_Populate(libzmq)
add_subdirectory(${libzmq_SOURCE_DIR} ${libzmq_BINARY_DIR} EXCLUDE_FROM_ALL)
endif()
else()
find_package(ZeroMQ 4 REQUIRED)
endif()
if (AARE_FETCH_FMT)
set(FMT_TEST OFF CACHE INTERNAL "disabling fmt tests")
FetchContent_Declare(
fmt
GIT_REPOSITORY https://github.com/fmtlib/fmt.git
GIT_TAG 10.2.1
GIT_PROGRESS TRUE
USES_TERMINAL_DOWNLOAD TRUE
)
set(FMT_INSTALL ON CACHE BOOL "")
# set(FMT_CMAKE_DIR "")
FetchContent_MakeAvailable(fmt)
set_property(TARGET fmt PROPERTY POSITION_INDEPENDENT_CODE ON)
install(TARGETS fmt
EXPORT ${project}-targets
LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
RUNTIME DESTINATION ${CMAKE_INSTALL_BINDIR}
INCLUDES DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}
)
else()
find_package(fmt 6 REQUIRED)
endif()
if (AARE_FETCH_JSON)
FetchContent_Declare(
json
URL https://github.com/nlohmann/json/releases/download/v3.11.3/json.tar.xz
)
set(JSON_Install ON CACHE BOOL "")
FetchContent_MakeAvailable(json)
set(NLOHMANN_JSON_TARGET_NAME nlohmann_json)
install(
TARGETS nlohmann_json
EXPORT "${TARGETS_EXPORT_NAME}"
)
message(STATUS "target: ${NLOHMANN_JSON_TARGET_NAME}")
else()
find_package(nlohmann_json 3.11.3 REQUIRED)
endif()
include(GNUInstallDirs)
# If conda build, always set lib dir to 'lib'
if($ENV{CONDA_BUILD})
set(CMAKE_INSTALL_LIBDIR "lib")
endif()
# Set lower / upper case project names
string(TOUPPER "${PROJECT_NAME}" PROJECT_NAME_UPPER)
string(TOLOWER "${PROJECT_NAME}" PROJECT_NAME_LOWER)
# Set targets export name (used by slsDetectorPackage and dependencies)
set(TARGETS_EXPORT_NAME "${PROJECT_NAME_LOWER}-targets")
set(namespace "aare::")
set(CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/cmake" ${CMAKE_MODULE_PATH})
# Check if project is being used directly or via add_subdirectory
set(AARE_MASTER_PROJECT OFF)
if (CMAKE_CURRENT_SOURCE_DIR STREQUAL CMAKE_SOURCE_DIR)
set(AARE_MASTER_PROJECT ON)
endif()
add_library(aare_compiler_flags INTERFACE)
target_compile_features(aare_compiler_flags INTERFACE cxx_std_17)
if(AARE_PYTHON_BINDINGS)
add_subdirectory(python)
endif()
#################
# MSVC specific #
#################
if(MSVC)
add_compile_definitions(AARE_MSVC)
if(CMAKE_BUILD_TYPE STREQUAL "Release")
message(STATUS "Release build")
target_compile_options(aare_compiler_flags INTERFACE /O2)
else()
message(STATUS "Debug build")
target_compile_options(
aare_compiler_flags
INTERFACE
/Od
/Zi
/MDd
/D_ITERATOR_DEBUG_LEVEL=2
)
target_link_options(
aare_compiler_flags
INTERFACE
/DEBUG:FULL
)
endif()
target_compile_options(
aare_compiler_flags
INTERFACE
/w # disable warnings
)
else()
######################
# GCC/Clang specific #
######################
if(CMAKE_BUILD_TYPE STREQUAL "Release")
message(STATUS "Release build")
target_compile_options(aare_compiler_flags INTERFACE -O3)
else()
message(STATUS "Debug build")
endif()
# Common flags for GCC and Clang
target_compile_options(
aare_compiler_flags
INTERFACE
-Wall
-Wextra
-pedantic
-Wshadow
-Wold-style-cast
-Wnon-virtual-dtor
-Woverloaded-virtual
-Wdouble-promotion
-Wformat=2
-Wredundant-decls
-Wvla
-Wdouble-promotion
-Werror=return-type #important can cause segfault in optimzed builds
)
endif() #GCC/Clang specific
if(AARE_ASAN)
message(STATUS "AddressSanitizer enabled")
target_compile_options(
aare_compiler_flags
INTERFACE
-fsanitize=address,undefined,pointer-compare
-fno-omit-frame-pointer
)
target_link_libraries(
aare_compiler_flags
INTERFACE
-fsanitize=address,undefined,pointer-compare
-fno-omit-frame-pointer
)
endif()
if(AARE_TESTS)
enable_testing()
add_subdirectory(tests)
endif()
###------------------------------------------------------------------------------MAIN LIBRARY
###------------------------------------------------------------------------------------------
set(PUBLICHEADERS
include/aare/ArrayExpr.hpp
include/aare/ClusterFinder.hpp
include/aare/ClusterFile.hpp
include/aare/CtbRawFile.hpp
include/aare/ClusterVector.hpp
include/aare/defs.hpp
include/aare/Dtype.hpp
include/aare/File.hpp
include/aare/FileInterface.hpp
include/aare/Frame.hpp
include/aare/NDArray.hpp
include/aare/NDView.hpp
include/aare/NumpyFile.hpp
include/aare/NumpyHelpers.hpp
include/aare/Pedestal.hpp
include/aare/PixelMap.hpp
include/aare/RawFile.hpp
include/aare/RawMasterFile.hpp
include/aare/RawSubFile.hpp
include/aare/VarClusterFinder.hpp
)
set(SourceFiles
${CMAKE_CURRENT_SOURCE_DIR}/src/CtbRawFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/defs.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/File.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/PixelMap.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/RawSubFile.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.cpp
)
add_library(aare_core STATIC ${SourceFiles})
target_include_directories(aare_core PUBLIC
"$<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>"
"$<INSTALL_INTERFACE:${CMAKE_INSTALL_INCLUDEDIR}>"
)
target_link_libraries(
aare_core
PUBLIC
fmt::fmt
nlohmann_json::nlohmann_json
${STD_FS_LIB} # from helpers.cmake
PRIVATE
aare_compiler_flags
)
set_target_properties(aare_core PROPERTIES
ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}
PUBLIC_HEADER "${PUBLICHEADERS}"
)
if (AARE_PYTHON_BINDINGS)
set_property(TARGET aare_core PROPERTY POSITION_INDEPENDENT_CODE ON)
endif()
if(AARE_TESTS)
set(TestSources
${CMAKE_CURRENT_SOURCE_DIR}/src/defs.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Dtype.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Frame.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/RawMasterFile.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NDArray.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NDView.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterFinder.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/ClusterVector.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/Pedestal.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyFile.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/NumpyHelpers.test.cpp
${CMAKE_CURRENT_SOURCE_DIR}/src/RawFile.test.cpp
)
target_sources(tests PRIVATE ${TestSources} )
endif()
###------------------------------------------------------------------------------------------
###------------------------------------------------------------------------------------------
if(AARE_MASTER_PROJECT)
install(TARGETS aare_core aare_compiler_flags
EXPORT "${TARGETS_EXPORT_NAME}"
LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
ARCHIVE DESTINATION ${CMAKE_INSTALL_LIBDIR}
PUBLIC_HEADER DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/aare
)
endif()
set(CMAKE_POSITION_INDEPENDENT_CODE ON)
set(CMAKE_INSTALL_RPATH $ORIGIN)
set(CMAKE_BUILD_WITH_INSTALL_RPATH FALSE)
# #Overall target to link to when using the library
# add_library(aare INTERFACE)
# target_link_libraries(aare INTERFACE aare_core aare_compiler_flags)
# target_include_directories(aare INTERFACE
# $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}/include>
# $<INSTALL_INTERFACE:include>
# )
# add_subdirectory(examples)
if(AARE_DOCS)
add_subdirectory(docs)
endif()
# custom target to run check formatting with clang-format
add_custom_target(
check-format
COMMAND find \( -name "*.cpp" -o -name "*.hpp" \) -not -path "./build/*" | xargs -I {} -n 1 -P 10 bash -c "clang-format -Werror -style=\"file:.clang-format\" {} | diff {} -"
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
COMMENT "Checking code formatting with clang-format"
VERBATIM
)
add_custom_target(
format-files
COMMAND find \( -name "*.cpp" -o -name "*.hpp" \) -not -path "./build/*" | xargs -I {} -n 1 -P 10 bash -c "clang-format -i -style=\"file:.clang-format\" {}"
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
COMMENT "Formatting with clang-format"
VERBATIM
)
if (AARE_IN_GITHUB_ACTIONS)
message(STATUS "Running in Github Actions")
set(CLANG_TIDY_COMMAND "clang-tidy-17")
else()
set(CLANG_TIDY_COMMAND "clang-tidy")
endif()
add_custom_target(
clang-tidy
COMMAND find \( -path "./src/*" -a -not -path "./src/python/*" -a \( -name "*.cpp" -not -name "*.test.cpp" \) \) -not -name "CircularFifo.hpp" -not -name "ProducerConsumerQueue.hpp" -not -name "VariableSizeClusterFinder.hpp" | xargs -I {} -n 1 -P 10 bash -c "${CLANG_TIDY_COMMAND} --config-file=.clang-tidy -p build {}"
WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
COMMENT "linting with clang-tidy"
VERBATIM
)
if(AARE_MASTER_PROJECT)
set(CMAKE_INSTALL_DIR "share/cmake/${PROJECT_NAME}")
set(PROJECT_LIBRARIES aare-core aare-compiler-flags )
include(cmake/package_config.cmake)
endif()

View File

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

11
benchmarks/CMakeLists.txt Normal file
View File

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

View File

@ -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();

11
cmake/FindSphinx.cmake Normal file
View File

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

46
cmake/helpers.cmake Normal file
View File

@ -0,0 +1,46 @@
function(default_build_type val)
if (NOT CMAKE_BUILD_TYPE AND NOT CMAKE_CONFIGURATION_TYPES)
message(STATUS "No build type selected, default to Release")
set(CMAKE_BUILD_TYPE ${val} CACHE STRING "Build type (default ${val})" FORCE)
endif()
endfunction()
function(set_std_fs_lib)
# from pybind11
# Check if we need to add -lstdc++fs or -lc++fs or nothing
if(DEFINED CMAKE_CXX_STANDARD AND CMAKE_CXX_STANDARD LESS 17)
set(STD_FS_NO_LIB_NEEDED TRUE)
elseif(MSVC)
set(STD_FS_NO_LIB_NEEDED TRUE)
else()
file(
WRITE ${CMAKE_CURRENT_BINARY_DIR}/main.cpp
"#include <filesystem>\nint main(int argc, char ** argv) {\n std::filesystem::path p(argv[0]);\n return p.string().length();\n}"
)
try_compile(
STD_FS_NO_LIB_NEEDED ${CMAKE_CURRENT_BINARY_DIR}
SOURCES ${CMAKE_CURRENT_BINARY_DIR}/main.cpp
COMPILE_DEFINITIONS -std=c++17)
try_compile(
STD_FS_NEEDS_STDCXXFS ${CMAKE_CURRENT_BINARY_DIR}
SOURCES ${CMAKE_CURRENT_BINARY_DIR}/main.cpp
COMPILE_DEFINITIONS -std=c++17
LINK_LIBRARIES stdc++fs)
try_compile(
STD_FS_NEEDS_CXXFS ${CMAKE_CURRENT_BINARY_DIR}
SOURCES ${CMAKE_CURRENT_BINARY_DIR}/main.cpp
COMPILE_DEFINITIONS -std=c++17
LINK_LIBRARIES c++fs)
endif()
if(${STD_FS_NEEDS_STDCXXFS})
set(STD_FS_LIB stdc++fs PARENT_SCOPE)
elseif(${STD_FS_NEEDS_CXXFS})
set(STD_FS_LIB c++fs PARENT_SCOPE)
elseif(${STD_FS_NO_LIB_NEEDED})
set(STD_FS_LIB "" PARENT_SCOPE)
else()
message(WARNING "Unknown C++17 compiler - not passing -lstdc++fs")
set(STD_FS_LIB "")
endif()
endfunction()

View File

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

View File

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

View File

@ -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

52
conda-recipe/meta.yaml Normal file
View File

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

77
docs/CMakeLists.txt Normal file
View File

@ -0,0 +1,77 @@
find_package(Doxygen REQUIRED)
find_package(Sphinx REQUIRED)
#Doxygen
set(DOXYGEN_IN ${CMAKE_CURRENT_SOURCE_DIR}/Doxyfile.in)
set(DOXYGEN_OUT ${CMAKE_CURRENT_BINARY_DIR}/Doxyfile)
configure_file(${DOXYGEN_IN} ${DOXYGEN_OUT} @ONLY)
#Sphinx
set(SPHINX_SOURCE ${CMAKE_CURRENT_SOURCE_DIR}/src)
set(SPHINX_BUILD ${CMAKE_CURRENT_BINARY_DIR})
file(GLOB SPHINX_SOURCE_FILES CONFIGURE_DEPENDS "src/*.rst")
# set(SPHINX_SOURCE_FILES
# src/index.rst
# src/Installation.rst
# src/Requirements.rst
# src/NDArray.rst
# src/NDView.rst
# src/File.rst
# src/Frame.rst
# src/Dtype.rst
# src/ClusterFinder.rst
# src/ClusterFile.rst
# src/Pedestal.rst
# src/RawFile.rst
# src/RawSubFile.rst
# src/RawMasterFile.rst
# src/VarClusterFinder.rst
# src/pyVarClusterFinder.rst
# src/pyFile.rst
# src/pyCtbRawFile.rst
# src/pyRawFile.rst
# src/pyRawMasterFile.rst
# )
foreach(filename ${SPHINX_SOURCE_FILES})
get_filename_component(fname ${filename} NAME)
message(STATUS "Copying ${filename} to ${SPHINX_BUILD}/src/${fname}")
configure_file(${filename} "${SPHINX_BUILD}/src/${fname}")
endforeach(filename ${SPHINX_SOURCE_FILES})
configure_file(
"${CMAKE_CURRENT_SOURCE_DIR}/conf.py.in"
"${SPHINX_BUILD}/conf.py"
@ONLY
)
configure_file(
"${CMAKE_CURRENT_SOURCE_DIR}/static/extra.css"
"${SPHINX_BUILD}/static/css/extra.css"
@ONLY
)
add_custom_target(
docs
COMMAND ${DOXYGEN_EXECUTABLE} ${DOXYGEN_OUT}
COMMAND ${SPHINX_EXECUTABLE} -a -b html
-Dbreathe_projects.aare=${CMAKE_CURRENT_BINARY_DIR}/xml
-c "${SPHINX_BUILD}"
${SPHINX_BUILD}/src
${SPHINX_BUILD}/html
COMMENT "Generating documentation with Sphinx"
)
add_custom_target(
rst
COMMAND ${SPHINX_EXECUTABLE} -a -b html
-Dbreathe_projects.aare=${CMAKE_CURRENT_BINARY_DIR}/xml
-c "${SPHINX_BUILD}"
${SPHINX_BUILD}/src
${SPHINX_BUILD}/html
COMMENT "Generating documentation with Sphinx"
)

1917
docs/Doxyfile.in Normal file

File diff suppressed because it is too large Load Diff

63
docs/conf.py.in Normal file
View File

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

7
docs/src/ClusterFile.rst Normal file
View File

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

View File

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

View File

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

19
docs/src/Consume.rst Normal file
View File

@ -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...

7
docs/src/Dtype.rst Normal file
View File

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

8
docs/src/File.rst Normal file
View File

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

8
docs/src/Frame.rst Normal file
View File

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

106
docs/src/Installation.rst Normal file
View File

@ -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.

7
docs/src/NDArray.rst Normal file
View File

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

7
docs/src/NDView.rst Normal file
View File

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

8
docs/src/Pedestal.rst Normal file
View File

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

8
docs/src/RawFile.rst Normal file
View File

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

View File

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

8
docs/src/RawSubFile.rst Normal file
View File

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

23
docs/src/Requirements.rst Normal file
View File

@ -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

View File

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

23
docs/src/_install.rst Normal file
View File

@ -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

58
docs/src/index.rst Normal file
View File

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

View File

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

11
docs/src/pyCtbRawFile.rst Normal file
View File

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

11
docs/src/pyFile.rst Normal file
View File

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

10
docs/src/pyRawFile.rst Normal file
View File

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

View File

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

View File

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

4
docs/static/extra.css vendored Normal file
View File

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

View File

@ -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

View File

@ -0,0 +1,81 @@
#pragma once
#include "aare/ClusterVector.hpp"
#include "aare/NDArray.hpp"
#include "aare/defs.hpp"
#include <filesystem>
#include <fstream>
namespace aare {
struct Cluster3x3 {
int16_t x;
int16_t y;
int32_t data[9];
};
typedef enum {
cBottomLeft = 0,
cBottomRight = 1,
cTopLeft = 2,
cTopRight = 3
} corner;
typedef enum {
pBottomLeft = 0,
pBottom = 1,
pBottomRight = 2,
pLeft = 3,
pCenter = 4,
pRight = 5,
pTopLeft = 6,
pTop = 7,
pTopRight = 8
} pixel;
struct ClusterAnalysis {
uint32_t c;
int32_t tot;
double etax;
double etay;
};
/*
Binary cluster file. Expects data to be layed out as:
int32_t frame_number
uint32_t number_of_clusters
int16_t x, int16_t y, int32_t data[9] x number_of_clusters
int32_t frame_number
uint32_t number_of_clusters
....
*/
class ClusterFile {
FILE *fp{};
uint32_t m_num_left{};
size_t m_chunk_size{};
const std::string m_mode;
public:
ClusterFile(const std::filesystem::path &fname, size_t chunk_size = 1000,
const std::string &mode = "r");
~ClusterFile();
std::vector<Cluster3x3> read_clusters(size_t n_clusters);
std::vector<Cluster3x3> read_frame(int32_t &out_fnum);
void write_frame(int32_t frame_number,
const ClusterVector<int32_t> &clusters);
std::vector<Cluster3x3>
read_cluster_with_cut(size_t n_clusters, double *noise_map, int nx, int ny);
size_t chunk_size() const { return m_chunk_size; }
void close();
};
int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad,
double *eta2x, double *eta2y, double *eta3x, double *eta3y);
int analyze_cluster(Cluster3x3& cl, int32_t *t2, int32_t *t3, char *quad,
double *eta2x, double *eta2y, double *eta3x, double *eta3y);
NDArray<double, 2> calculate_eta2( ClusterVector<int>& clusters);
std::array<double,2> calculate_eta2( Cluster3x3& cl);
} // namespace aare

View File

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

View File

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

View File

@ -0,0 +1,180 @@
#pragma once
#include <cstddef>
#include <cstdint>
#include <numeric>
#include <vector>
#include <fmt/core.h>
namespace aare {
/**
* @brief ClusterVector is a container for clusters of various sizes. It uses a
* contiguous memory buffer to store the clusters.
* @note push_back can invalidate pointers to elements in the container
* @tparam T data type of the pixels in the cluster
* @tparam CoordType data type of the x and y coordinates of the cluster (normally int16_t)
*/
template <typename T, typename CoordType=int16_t> class ClusterVector {
using value_type = T;
size_t m_cluster_size_x;
size_t m_cluster_size_y;
std::byte *m_data{};
size_t m_size{0};
size_t m_capacity;
/*
Format string used in the python bindings to create a numpy
array from the buffer
= - native byte order
h - short
d - double
i - int
*/
constexpr static char m_fmt_base[] = "=h:x:\nh:y:\n({},{}){}:data:" ;
public:
/**
* @brief Construct a new ClusterVector object
* @param cluster_size_x size of the cluster in x direction
* @param cluster_size_y size of the cluster in y direction
* @param capacity initial capacity of the buffer in number of clusters
*/
ClusterVector(size_t cluster_size_x, size_t cluster_size_y,
size_t capacity = 1024)
: m_cluster_size_x(cluster_size_x), m_cluster_size_y(cluster_size_y),
m_capacity(capacity) {
allocate_buffer(capacity);
}
~ClusterVector() {
delete[] m_data;
}
//Move constructor
ClusterVector(ClusterVector &&other) noexcept
: m_cluster_size_x(other.m_cluster_size_x),
m_cluster_size_y(other.m_cluster_size_y), m_data(other.m_data),
m_size(other.m_size), m_capacity(other.m_capacity) {
other.m_data = nullptr;
other.m_size = 0;
other.m_capacity = 0;
}
//Move assignment operator
ClusterVector& operator=(ClusterVector &&other) noexcept {
if (this != &other) {
delete[] m_data;
m_cluster_size_x = other.m_cluster_size_x;
m_cluster_size_y = other.m_cluster_size_y;
m_data = other.m_data;
m_size = other.m_size;
m_capacity = other.m_capacity;
other.m_data = nullptr;
other.m_size = 0;
other.m_capacity = 0;
}
return *this;
}
/**
* @brief Reserve space for at least capacity clusters
* @param capacity number of clusters to reserve space for
* @note If capacity is less than the current capacity, the function does nothing.
*/
void reserve(size_t capacity) {
if (capacity > m_capacity) {
allocate_buffer(capacity);
}
}
/**
* @brief Add a cluster to the vector
* @param x x-coordinate of the cluster
* @param y y-coordinate of the cluster
* @param data pointer to the data of the cluster
* @warning The data pointer must point to a buffer of size cluster_size_x * cluster_size_y * sizeof(T)
*/
void push_back(CoordType x, CoordType y, const std::byte *data) {
if (m_size == m_capacity) {
allocate_buffer(m_capacity * 2);
}
std::byte *ptr = element_ptr(m_size);
*reinterpret_cast<CoordType *>(ptr) = x;
ptr += sizeof(CoordType);
*reinterpret_cast<CoordType *>(ptr) = y;
ptr += sizeof(CoordType);
std::copy(data, data + m_cluster_size_x * m_cluster_size_y * sizeof(T),
ptr);
m_size++;
}
/**
* @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 = element_offset();
const size_t n_pixels = m_cluster_size_x * m_cluster_size_y;
std::byte *ptr = m_data + 2 * sizeof(CoordType); // skip x and y
for (size_t i = 0; i < m_size; i++) {
sums[i] =
std::accumulate(reinterpret_cast<T *>(ptr),
reinterpret_cast<T *>(ptr) + n_pixels, T{});
ptr += stride;
}
return sums;
}
size_t size() const { return m_size; }
size_t capacity() const { return m_capacity; }
/**
* @brief Return the offset in bytes for a single cluster
*/
size_t element_offset() const {
return 2*sizeof(CoordType) +
m_cluster_size_x * m_cluster_size_y * sizeof(T);
}
/**
* @brief Return the offset in bytes for the i-th cluster
*/
size_t element_offset(size_t i) const { return element_offset() * i; }
/**
* @brief Return a pointer to the i-th cluster
*/
std::byte *element_ptr(size_t i) { return m_data + element_offset(i); }
const std::byte * element_ptr(size_t i) const { return m_data + element_offset(i); }
size_t cluster_size_x() const { return m_cluster_size_x; }
size_t cluster_size_y() const { return m_cluster_size_y; }
std::byte *data() { return m_data; }
std::byte const *data() const { return m_data; }
template<typename V>
V& at(size_t i) {
return *reinterpret_cast<V*>(element_ptr(i));
}
const std::string_view fmt_base() const {
//TODO! how do we match on coord_t?
return m_fmt_base;
}
private:
void allocate_buffer(size_t new_capacity) {
size_t num_bytes = element_offset() * new_capacity;
std::byte *new_data = new std::byte[num_bytes]{};
std::copy(m_data, m_data + element_offset() * m_size, new_data);
delete[] m_data;
m_data = new_data;
m_capacity = new_capacity;
}
};
} // namespace aare

View File

@ -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);
};
}

83
include/aare/Dtype.hpp Normal file
View File

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

64
include/aare/File.hpp Normal file
View File

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

View File

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

124
include/aare/Frame.hpp Normal file
View File

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

421
include/aare/NDArray.hpp Normal file
View File

@ -0,0 +1,421 @@
#pragma once
/*
Container holding image data, or a time series of image data in contigious
memory.
TODO! Add expression templates for operators
*/
#include "aare/ArrayExpr.hpp"
#include "aare/NDView.hpp"
#include <algorithm>
#include <array>
#include <cmath>
#include <fmt/format.h>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <numeric>
namespace aare {
template <typename T, int64_t Ndim = 2>
class NDArray : public ArrayExpr<NDArray<T, Ndim>, Ndim> {
std::array<int64_t, Ndim> shape_;
std::array<int64_t, Ndim> strides_;
size_t size_{};
T *data_;
public:
/**
* @brief Default constructor. Will construct an empty NDArray.
*
*/
NDArray() : shape_(), strides_(c_strides<Ndim>(shape_)), data_(nullptr) {};
/**
* @brief Construct a new NDArray object with a given shape.
* @note The data is uninitialized.
*
* @param shape shape of the new NDArray
*/
explicit NDArray(std::array<int64_t, Ndim> shape)
: shape_(shape), strides_(c_strides<Ndim>(shape_)),
size_(std::accumulate(shape_.begin(), shape_.end(), 1,
std::multiplies<>())),
data_(new T[size_]) {}
/**
* @brief Construct a new NDArray object with a shape and value.
*
* @param shape shape of the new array
* @param value value to initialize the array with
*/
NDArray(std::array<int64_t, Ndim> shape, T value) : NDArray(shape) {
this->operator=(value);
}
/**
* @brief Construct a new NDArray object from a NDView.
* @note The data is copied from the view to the NDArray.
*
* @param v view of data to initialize the NDArray with
*/
explicit NDArray(const NDView<T, Ndim> v) : NDArray(v.shape()) {
std::copy(v.begin(), v.end(), begin());
}
// Move constructor
NDArray(NDArray &&other) noexcept
: shape_(other.shape_), strides_(c_strides<Ndim>(shape_)),
size_(other.size_), data_(other.data_) {
other.reset(); // TODO! is this necessary?
}
// Copy constructor
NDArray(const NDArray &other)
: shape_(other.shape_), strides_(c_strides<Ndim>(shape_)),
size_(other.size_), data_(new T[size_]) {
std::copy(other.data_, other.data_ + size_, data_);
}
// Conversion operator from array expression to array
template <typename E>
NDArray(ArrayExpr<E, Ndim> &&expr) : NDArray(expr.shape()) {
for (size_t i = 0; i < size_; ++i) {
data_[i] = expr[i];
}
}
~NDArray() { delete[] data_; }
auto begin() { return data_; }
auto end() { return data_ + size_; }
using value_type = T;
NDArray &operator=(NDArray &&other) noexcept; // Move assign
NDArray &operator=(const NDArray &other); // Copy assign
NDArray &operator+=(const NDArray &other);
NDArray &operator-=(const NDArray &other);
NDArray &operator*=(const NDArray &other);
// NDArray& operator/=(const NDArray& other);
template <typename V> NDArray &operator/=(const NDArray<V, Ndim> &other) {
// check shape
if (shape_ == other.shape()) {
for (uint32_t i = 0; i < size_; ++i) {
data_[i] /= other(i);
}
return *this;
}
throw(std::runtime_error("Shape of NDArray must match"));
}
NDArray<bool, Ndim> operator>(const NDArray &other);
bool operator==(const NDArray &other) const;
bool operator!=(const NDArray &other) const;
NDArray &operator=(const T & /*value*/);
NDArray &operator+=(const T & /*value*/);
NDArray operator+(const T & /*value*/);
NDArray &operator-=(const T & /*value*/);
NDArray operator-(const T & /*value*/);
NDArray &operator*=(const T & /*value*/);
NDArray operator*(const T & /*value*/);
NDArray &operator/=(const T & /*value*/);
NDArray operator/(const T & /*value*/);
NDArray &operator&=(const T & /*mask*/);
void sqrt() {
for (int i = 0; i < size_; ++i) {
data_[i] = std::sqrt(data_[i]);
}
}
NDArray &operator++(); // pre inc
template <typename... Ix>
std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) {
return data_[element_offset(strides_, index...)];
}
template <typename... Ix>
std::enable_if_t<sizeof...(Ix) == Ndim, T &> operator()(Ix... index) const {
return data_[element_offset(strides_, index...)];
}
template <typename... Ix>
std::enable_if_t<sizeof...(Ix) == Ndim, T> value(Ix... index) {
return data_[element_offset(strides_, index...)];
}
// TODO! is int the right type for index?
T &operator()(int64_t i) { return data_[i]; }
const T &operator()(int64_t i) const { return data_[i]; }
T &operator[](int64_t i) { return data_[i]; }
const T &operator[](int64_t i) const { return data_[i]; }
T *data() { return data_; }
std::byte *buffer() { return reinterpret_cast<std::byte *>(data_); }
size_t size() const { return size_; }
size_t total_bytes() const { return size_ * sizeof(T); }
std::array<int64_t, Ndim> shape() const noexcept { return shape_; }
int64_t shape(int64_t i) const noexcept { return shape_[i]; }
std::array<int64_t, Ndim> strides() const noexcept { return strides_; }
size_t bitdepth() const noexcept { return sizeof(T) * 8; }
std::array<int64_t, Ndim> byte_strides() const noexcept {
auto byte_strides = strides_;
for (auto &val : byte_strides)
val *= sizeof(T);
return byte_strides;
}
/**
* @brief Create a view of the NDArray.
*
* @return NDView<T, Ndim>
*/
NDView<T, Ndim> view() const { return NDView<T, Ndim>{data_, shape_}; }
void Print();
void Print_all();
void Print_some();
void reset() {
data_ = nullptr;
size_ = 0;
std::fill(shape_.begin(), shape_.end(), 0);
std::fill(strides_.begin(), strides_.end(), 0);
}
};
// Move assign
template <typename T, int64_t Ndim>
NDArray<T, Ndim> &
NDArray<T, Ndim>::operator=(NDArray<T, Ndim> &&other) noexcept {
if (this != &other) {
delete[] data_;
data_ = other.data_;
shape_ = other.shape_;
size_ = other.size_;
strides_ = other.strides_;
other.reset();
}
return *this;
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const NDArray<T, Ndim> &other) {
// check shape
if (shape_ == other.shape_) {
for (size_t i = 0; i < size_; ++i) {
data_[i] += other.data_[i];
}
return *this;
}
throw(std::runtime_error("Shape of ImageDatas must match"));
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator-=(const NDArray<T, Ndim> &other) {
// check shape
if (shape_ == other.shape_) {
for (uint32_t i = 0; i < size_; ++i) {
data_[i] -= other.data_[i];
}
return *this;
}
throw(std::runtime_error("Shape of ImageDatas must match"));
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator*=(const NDArray<T, Ndim> &other) {
// check shape
if (shape_ == other.shape_) {
for (uint32_t i = 0; i < size_; ++i) {
data_[i] *= other.data_[i];
}
return *this;
}
throw(std::runtime_error("Shape of ImageDatas must match"));
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator&=(const T &mask) {
for (auto it = begin(); it != end(); ++it)
*it &= mask;
return *this;
}
template <typename T, int64_t Ndim>
NDArray<bool, Ndim> NDArray<T, Ndim>::operator>(const NDArray &other) {
if (shape_ == other.shape_) {
NDArray<bool, Ndim> result{shape_};
for (int i = 0; i < size_; ++i) {
result(i) = (data_[i] > other.data_[i]);
}
return result;
}
throw(std::runtime_error("Shape of ImageDatas must match"));
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator=(const NDArray<T, Ndim> &other) {
if (this != &other) {
delete[] data_;
shape_ = other.shape_;
strides_ = other.strides_;
size_ = other.size_;
data_ = new T[size_];
std::copy(other.data_, other.data_ + size_, data_);
}
return *this;
}
template <typename T, int64_t Ndim>
bool NDArray<T, Ndim>::operator==(const NDArray<T, Ndim> &other) const {
if (shape_ != other.shape_)
return false;
for (uint32_t i = 0; i != size_; ++i)
if (data_[i] != other.data_[i])
return false;
return true;
}
template <typename T, int64_t Ndim>
bool NDArray<T, Ndim>::operator!=(const NDArray<T, Ndim> &other) const {
return !((*this) == other);
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator++() {
for (uint32_t i = 0; i < size_; ++i)
data_[i] += 1;
return *this;
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator=(const T &value) {
std::fill_n(data_, size_, value);
return *this;
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator+=(const T &value) {
for (uint32_t i = 0; i < size_; ++i)
data_[i] += value;
return *this;
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator+(const T &value) {
NDArray result = *this;
result += value;
return result;
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator-=(const T &value) {
for (uint32_t i = 0; i < size_; ++i)
data_[i] -= value;
return *this;
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator-(const T &value) {
NDArray result = *this;
result -= value;
return result;
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator/=(const T &value) {
for (uint32_t i = 0; i < size_; ++i)
data_[i] /= value;
return *this;
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator/(const T &value) {
NDArray result = *this;
result /= value;
return result;
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> &NDArray<T, Ndim>::operator*=(const T &value) {
for (uint32_t i = 0; i < size_; ++i)
data_[i] *= value;
return *this;
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> NDArray<T, Ndim>::operator*(const T &value) {
NDArray result = *this;
result *= value;
return result;
}
template <typename T, int64_t Ndim> void NDArray<T, Ndim>::Print() {
if (shape_[0] < 20 && shape_[1] < 20)
Print_all();
else
Print_some();
}
template <typename T, int64_t Ndim>
std::ostream &operator<<(std::ostream &os, const NDArray<T, Ndim> &arr) {
for (auto row = 0; row < arr.shape(0); ++row) {
for (auto col = 0; col < arr.shape(1); ++col) {
os << std::setw(3);
os << arr(row, col) << " ";
}
os << "\n";
}
return os;
}
template <typename T, int64_t Ndim> void NDArray<T, Ndim>::Print_all() {
for (auto row = 0; row < shape_[0]; ++row) {
for (auto col = 0; col < shape_[1]; ++col) {
std::cout << std::setw(3);
std::cout << (*this)(row, col) << " ";
}
std::cout << "\n";
}
}
template <typename T, int64_t Ndim> void NDArray<T, Ndim>::Print_some() {
for (auto row = 0; row < 5; ++row) {
for (auto col = 0; col < 5; ++col) {
std::cout << std::setw(7);
std::cout << (*this)(row, col) << " ";
}
std::cout << "\n";
}
}
template <typename T, int64_t Ndim>
void save(NDArray<T, Ndim> &img, std::string &pathname) {
std::ofstream f;
f.open(pathname, std::ios::binary);
f.write(img.buffer(), img.size() * sizeof(T));
f.close();
}
template <typename T, int64_t Ndim>
NDArray<T, Ndim> load(const std::string &pathname,
std::array<int64_t, Ndim> shape) {
NDArray<T, Ndim> img{shape};
std::ifstream f;
f.open(pathname, std::ios::binary);
f.read(img.buffer(), img.size() * sizeof(T));
f.close();
return img;
}
} // namespace aare

178
include/aare/NDView.hpp Normal file
View File

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

119
include/aare/NumpyFile.hpp Normal file
View File

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

View File

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

207
include/aare/Pedestal.hpp Normal file
View File

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

20
include/aare/PixelMap.hpp Normal file
View File

@ -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

121
include/aare/RawFile.hpp Normal file
View File

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

View File

@ -0,0 +1,153 @@
#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();
};
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; }
};
/**
* @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

View File

@ -0,0 +1,72 @@
#pragma once
#include "aare/Frame.hpp"
#include "aare/defs.hpp"
#include <cstdint>
#include <filesystem>
#include <map>
#include <optional>
namespace aare {
/**
* @brief Class to read a singe subfile written in .raw format. Used from RawFile to read
* the entire detector. Can be used directly to read part of the image.
*/
class RawSubFile {
protected:
std::ifstream m_file;
DetectorType m_detector_type;
size_t m_bitdepth;
std::filesystem::path m_fname;
size_t m_rows{};
size_t m_cols{};
size_t m_bytes_per_frame{};
size_t n_frames{};
uint32_t m_pos_row{};
uint32_t m_pos_col{};
std::optional<NDArray<ssize_t, 2>> m_pixel_map;
public:
/**
* @brief SubFile constructor
* @param fname path to the subfile
* @param detector detector type
* @param rows number of rows in the subfile
* @param cols number of columns in the subfile
* @param bitdepth bitdepth of the subfile
* @throws std::invalid_argument if the detector,type pair is not supported
*/
RawSubFile(const std::filesystem::path &fname, DetectorType detector,
size_t rows, size_t cols, size_t bitdepth, uint32_t pos_row = 0, uint32_t pos_col = 0);
~RawSubFile() = default;
/**
* @brief Seek to the given frame number
* @note Puts the file pointer at the start of the header, not the start of the data
* @param frame_index frame position in file to seek to
* @throws std::runtime_error if the frame number is out of range
*/
void seek(size_t frame_index);
size_t tell();
void read_into(std::byte *image_buf, DetectorHeader *header = nullptr);
void get_part(std::byte *buffer, size_t frame_index);
void read_header(DetectorHeader *header);
size_t rows() const;
size_t cols() const;
size_t frame_number(size_t frame_index);
size_t bytes_per_frame() const { return m_bytes_per_frame; }
size_t pixels_per_frame() const { return m_rows * m_cols; }
size_t bytes_per_pixel() const { return m_bitdepth / 8; }
};
} // namespace aare

View File

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

233
include/aare/defs.hpp Normal file
View File

@ -0,0 +1,233 @@
#pragma once
#include "aare/Dtype.hpp"
// #include "aare/utils/logger.hpp"
#include <array>
#include <stdexcept>
#include <cassert>
#include <cstdint>
#include <cstring>
#include <string>
#include <string_view>
#include <variant>
#include <vector>
/**
* @brief LOCATION macro to get the current location in the code
*/
#define LOCATION \
std::string(__FILE__) + std::string(":") + std::to_string(__LINE__) + \
":" + std::string(__func__) + ":"
#ifdef AARE_CUSTOM_ASSERT
#define AARE_ASSERT(expr)\
if (expr)\
{}\
else\
aare::assert_failed(LOCATION + " Assertion failed: " + #expr + "\n");
#else
#define AARE_ASSERT(cond)\
do { (void)sizeof(cond); } while(0)
#endif
namespace aare {
void assert_failed(const std::string &msg);
class DynamicCluster {
public:
int cluster_sizeX;
int cluster_sizeY;
int16_t x;
int16_t y;
Dtype dt; // 4 bytes
private:
std::byte *m_data;
public:
DynamicCluster(int cluster_sizeX_, int cluster_sizeY_,
Dtype dt_ = Dtype(typeid(int32_t)))
: cluster_sizeX(cluster_sizeX_), cluster_sizeY(cluster_sizeY_),
dt(dt_) {
m_data = new std::byte[cluster_sizeX * cluster_sizeY * dt.bytes()]{};
}
DynamicCluster() : DynamicCluster(3, 3) {}
DynamicCluster(const DynamicCluster &other)
: DynamicCluster(other.cluster_sizeX, other.cluster_sizeY, other.dt) {
if (this == &other)
return;
x = other.x;
y = other.y;
memcpy(m_data, other.m_data, other.bytes());
}
DynamicCluster &operator=(const DynamicCluster &other) {
if (this == &other)
return *this;
this->~DynamicCluster();
new (this) DynamicCluster(other);
return *this;
}
DynamicCluster(DynamicCluster &&other) noexcept
: cluster_sizeX(other.cluster_sizeX),
cluster_sizeY(other.cluster_sizeY), x(other.x), y(other.y),
dt(other.dt), m_data(other.m_data) {
other.m_data = nullptr;
other.dt = Dtype(Dtype::TypeIndex::ERROR);
}
~DynamicCluster() { delete[] m_data; }
template <typename T> T get(int idx) {
(sizeof(T) == dt.bytes())
? 0
: throw std::invalid_argument("[ERROR] Type size mismatch");
return *reinterpret_cast<T *>(m_data + idx * dt.bytes());
}
template <typename T> auto set(int idx, T val) {
(sizeof(T) == dt.bytes())
? 0
: throw std::invalid_argument("[ERROR] Type size mismatch");
return memcpy(m_data + idx * dt.bytes(), &val, dt.bytes());
}
template <typename T> std::string to_string() const {
(sizeof(T) == dt.bytes())
? 0
: throw std::invalid_argument("[ERROR] Type size mismatch");
std::string s = "x: " + std::to_string(x) + " y: " + std::to_string(y) +
"\nm_data: [";
for (int i = 0; i < cluster_sizeX * cluster_sizeY; i++) {
s += std::to_string(
*reinterpret_cast<T *>(m_data + i * dt.bytes())) +
" ";
}
s += "]";
return s;
}
/**
* @brief size of the cluster in bytes when saved to a file
*/
size_t size() const { return cluster_sizeX * cluster_sizeY; }
size_t bytes() const { return cluster_sizeX * cluster_sizeY * dt.bytes(); }
auto begin() const { return m_data; }
auto end() const {
return m_data + cluster_sizeX * cluster_sizeY * dt.bytes();
}
std::byte *data() { return m_data; }
};
/**
* @brief header contained in parts of frames
*/
struct DetectorHeader {
uint64_t frameNumber;
uint32_t expLength;
uint32_t packetNumber;
uint64_t bunchId;
uint64_t timestamp;
uint16_t modId;
uint16_t row;
uint16_t column;
uint16_t reserved;
uint32_t debug;
uint16_t roundRNumber;
uint8_t detType;
uint8_t version;
std::array<uint8_t, 64> packetMask;
std::string to_string() {
std::string packetMaskStr = "[";
for (auto &i : packetMask) {
packetMaskStr += std::to_string(i) + ", ";
}
packetMaskStr += "]";
return "frameNumber: " + std::to_string(frameNumber) + "\n" +
"expLength: " + std::to_string(expLength) + "\n" +
"packetNumber: " + std::to_string(packetNumber) + "\n" +
"bunchId: " + std::to_string(bunchId) + "\n" +
"timestamp: " + std::to_string(timestamp) + "\n" +
"modId: " + std::to_string(modId) + "\n" +
"row: " + std::to_string(row) + "\n" +
"column: " + std::to_string(column) + "\n" +
"reserved: " + std::to_string(reserved) + "\n" +
"debug: " + std::to_string(debug) + "\n" +
"roundRNumber: " + std::to_string(roundRNumber) + "\n" +
"detType: " + std::to_string(detType) + "\n" +
"version: " + std::to_string(version) + "\n" +
"packetMask: " + packetMaskStr + "\n";
}
};
template <typename T> struct t_xy {
T row;
T col;
bool operator==(const t_xy &other) const {
return row == other.row && col == other.col;
}
bool operator!=(const t_xy &other) const { return !(*this == other); }
std::string to_string() const {
return "{ x: " + std::to_string(row) + " y: " + std::to_string(col) +
" }";
}
};
using xy = t_xy<uint32_t>;
struct ModuleGeometry{
int x{};
int y{};
int height{};
int width{};
};
using dynamic_shape = std::vector<int64_t>;
//TODO! Can we uniform enums between the libraries?
/**
* @brief Enum class to identify different detectors.
* The values are the same as in slsDetectorPackage
* Different spelling to avoid confusion with the slsDetectorPackage
*/
enum class DetectorType {
//Standard detectors match the enum values from slsDetectorPackage
Generic,
Eiger,
Gotthard,
Jungfrau,
ChipTestBoard,
Moench,
Mythen3,
Gotthard2,
Xilinx_ChipTestBoard,
//Additional detectors used for defining processing. Variants of the standard ones.
Moench03=100,
Moench03_old,
Unknown
};
enum class TimingMode { Auto, Trigger };
enum class FrameDiscardPolicy { NoDiscard, Discard, DiscardPartial };
template <class T> T StringTo(const std::string &arg) { return T(arg); }
template <class T> std::string ToString(T arg) { return T(arg); }
template <> DetectorType StringTo(const std::string & /*name*/);
template <> std::string ToString(DetectorType arg);
template <> TimingMode StringTo(const std::string & /*mode*/);
template <> FrameDiscardPolicy StringTo(const std::string & /*mode*/);
using DataTypeVariants = std::variant<uint16_t, uint32_t>;
} // namespace aare

15
pyproject.toml Normal file
View File

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

59
python/CMakeLists.txt Normal file
View File

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

191
python/aare/CtbRawFile.py Normal file
View File

@ -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

66
python/aare/RawFile.py Normal file
View File

@ -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

View File

@ -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__()

15
python/aare/__init__.py Normal file
View File

@ -0,0 +1,15 @@
# Make the compiled classes that live in _aare available from aare.
from . import _aare
from ._aare import File, RawMasterFile, RawSubFile
from ._aare import Pedestal_d, Pedestal_f, ClusterFinder, VarClusterFinder
from ._aare import DetectorType
from ._aare import ClusterFile
from ._aare import hitmap
from .CtbRawFile import CtbRawFile
from .RawFile import RawFile
from .ScanParameters import ScanParameters
from .utils import random_pixels, random_pixel

48
python/aare/transform.py Normal file
View File

@ -0,0 +1,48 @@
import numpy as np
from . import _aare
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()

23
python/aare/utils.py Normal file
View File

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

58
python/examples/play.py Normal file
View File

@ -0,0 +1,58 @@
import sys
sys.path.append('/home/l_msdetect/erik/aare/build')
#Our normal python imports
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
import boost_histogram as bh
import time
from aare import File, ClusterFinder, VarClusterFinder
base = Path('/mnt/sls_det_storage/matterhorn_data/aare_test_data/')
f = File(base/'Moench03new/cu_half_speed_master_4.json')
cf = ClusterFinder((400,400), (3,3))
for i in range(1000):
cf.push_pedestal_frame(f.read_frame())
fig, ax = plt.subplots()
im = ax.imshow(cf.pedestal())
cf.pedestal()
cf.noise()
N = 500
t0 = time.perf_counter()
hist1 = bh.Histogram(bh.axis.Regular(40, -2, 4000))
f.seek(0)
t0 = time.perf_counter()
data = f.read_n(N)
t_elapsed = time.perf_counter()-t0
n_bytes = data.itemsize*data.size
print(f'Reading {N} frames took {t_elapsed:.3f}s {N/t_elapsed:.0f} FPS, {n_bytes/1024**2:.4f} GB/s')
for frame in data:
a = cf.find_clusters(frame)
clusters = cf.steal_clusters()
# t_elapsed = time.perf_counter()-t0
# print(f'Clustering {N} frames took {t_elapsed:.2f}s {N/t_elapsed:.0f} FPS')
# t0 = time.perf_counter()
# total_clusters = clusters.size
# hist1.fill(clusters.sum())
# t_elapsed = time.perf_counter()-t0
# print(f'Filling histogram with the sum of {total_clusters} clusters took: {t_elapsed:.3f}s, {total_clusters/t_elapsed:.3g} clust/s')
# print(f'Average number of clusters per frame {total_clusters/N:.3f}')

124
python/src/cluster.hpp Normal file
View File

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

View File

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

View File

@ -0,0 +1,57 @@
#include "aare/CtbRawFile.hpp"
#include "aare/File.hpp"
#include "aare/Frame.hpp"
#include "aare/RawFile.hpp"
#include "aare/RawMasterFile.hpp"
#include "aare/RawSubFile.hpp"
#include "aare/defs.hpp"
// #include "aare/fClusterFileV2.hpp"
#include <cstdint>
#include <filesystem>
#include <pybind11/iostream.h>
#include <pybind11/numpy.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/stl/filesystem.h>
#include <string>
namespace py = pybind11;
using namespace ::aare;
void define_ctb_raw_file_io_bindings(py::module &m) {
py::class_<CtbRawFile>(m, "CtbRawFile")
.def(py::init<const std::filesystem::path &>())
.def("read_frame",
[](CtbRawFile &self) {
size_t image_size = self.image_size_in_bytes();
py::array image;
std::vector<ssize_t> shape;
shape.reserve(2);
shape.push_back(1);
shape.push_back(image_size);
py::array_t<DetectorHeader> header(1);
// always read bytes
image = py::array_t<uint8_t>(shape);
self.read_into(
reinterpret_cast<std::byte *>(image.mutable_data()),
header.mutable_data());
return py::make_tuple(header, image);
})
.def("seek", &CtbRawFile::seek)
.def("tell", &CtbRawFile::tell)
.def("master", &CtbRawFile::master)
.def_property_readonly("image_size_in_bytes",
&CtbRawFile::image_size_in_bytes)
.def_property_readonly("frames_in_file", &CtbRawFile::frames_in_file);
}

257
python/src/file.hpp Normal file
View File

@ -0,0 +1,257 @@
#include "aare/CtbRawFile.hpp"
#include "aare/File.hpp"
#include "aare/Frame.hpp"
#include "aare/RawFile.hpp"
#include "aare/RawMasterFile.hpp"
#include "aare/RawSubFile.hpp"
#include "aare/defs.hpp"
// #include "aare/fClusterFileV2.hpp"
#include <cstdint>
#include <filesystem>
#include <pybind11/iostream.h>
#include <pybind11/numpy.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
#include <pybind11/stl/filesystem.h>
#include <string>
namespace py = pybind11;
using namespace ::aare;
void define_file_io_bindings(py::module &m) {
py::enum_<DetectorType>(m, "DetectorType")
.value("Jungfrau", DetectorType::Jungfrau)
.value("Eiger", DetectorType::Eiger)
.value("Mythen3", DetectorType::Mythen3)
.value("Moench", DetectorType::Moench)
.value("Moench03", DetectorType::Moench03)
.value("Moench03_old", DetectorType::Moench03_old)
.value("ChipTestBoard", DetectorType::ChipTestBoard)
.value("Unknown", DetectorType::Unknown);
PYBIND11_NUMPY_DTYPE(DetectorHeader, frameNumber, expLength, packetNumber,
bunchId, timestamp, modId, row, column, reserved,
debug, roundRNumber, detType, version, packetMask);
py::class_<File>(m, "File")
.def(py::init([](const std::filesystem::path &fname) {
return File(fname, "r", {});
}))
.def(py::init(
[](const std::filesystem::path &fname, const std::string &mode) {
return File(fname, mode, {});
}))
.def(py::init<const std::filesystem::path &, const std::string &,
const FileConfig &>())
.def("frame_number", py::overload_cast<>(&File::frame_number))
.def("frame_number", py::overload_cast<size_t>(&File::frame_number))
.def_property_readonly("bytes_per_frame", &File::bytes_per_frame)
.def_property_readonly("pixels_per_frame", &File::pixels_per_frame)
.def("seek", &File::seek)
.def("tell", &File::tell)
.def_property_readonly("total_frames", &File::total_frames)
.def_property_readonly("rows", &File::rows)
.def_property_readonly("cols", &File::cols)
.def_property_readonly("bitdepth", &File::bitdepth)
.def_property_readonly("bytes_per_pixel", &File::bytes_per_pixel)
.def_property_readonly(
"detector_type",
[](File &self) { return ToString(self.detector_type()); })
.def("read_frame",
[](File &self) {
const uint8_t item_size = self.bytes_per_pixel();
py::array image;
std::vector<ssize_t> shape;
shape.reserve(2);
shape.push_back(self.rows());
shape.push_back(self.cols());
if (item_size == 1) {
image = py::array_t<uint8_t>(shape);
} else if (item_size == 2) {
image = py::array_t<uint16_t>(shape);
} else if (item_size == 4) {
image = py::array_t<uint32_t>(shape);
}
self.read_into(
reinterpret_cast<std::byte *>(image.mutable_data()));
return image;
})
.def("read_frame",
[](File &self, size_t frame_number) {
self.seek(frame_number);
const uint8_t item_size = self.bytes_per_pixel();
py::array image;
std::vector<ssize_t> shape;
shape.reserve(2);
shape.push_back(self.rows());
shape.push_back(self.cols());
if (item_size == 1) {
image = py::array_t<uint8_t>(shape);
} else if (item_size == 2) {
image = py::array_t<uint16_t>(shape);
} else if (item_size == 4) {
image = py::array_t<uint32_t>(shape);
}
self.read_into(
reinterpret_cast<std::byte *>(image.mutable_data()));
return image;
})
.def("read_n", [](File &self, size_t n_frames) {
//adjust for actual frames left in the file
n_frames = std::min(n_frames, self.total_frames()-self.tell());
if(n_frames == 0){
throw std::runtime_error("No frames left in file");
}
std::vector<size_t> shape{n_frames, self.rows(), self.cols()};
py::array image;
const uint8_t item_size = self.bytes_per_pixel();
if (item_size == 1) {
image = py::array_t<uint8_t>(shape);
} else if (item_size == 2) {
image = py::array_t<uint16_t>(shape);
} else if (item_size == 4) {
image = py::array_t<uint32_t>(shape);
}
self.read_into(reinterpret_cast<std::byte *>(image.mutable_data()),
n_frames);
return image;
});
py::class_<FileConfig>(m, "FileConfig")
.def(py::init<>())
.def_readwrite("rows", &FileConfig::rows)
.def_readwrite("cols", &FileConfig::cols)
.def_readwrite("version", &FileConfig::version)
.def_readwrite("geometry", &FileConfig::geometry)
.def_readwrite("detector_type", &FileConfig::detector_type)
.def_readwrite("max_frames_per_file", &FileConfig::max_frames_per_file)
.def_readwrite("total_frames", &FileConfig::total_frames)
.def_readwrite("dtype", &FileConfig::dtype)
.def("__eq__", &FileConfig::operator==)
.def("__ne__", &FileConfig::operator!=)
.def("__repr__", [](const FileConfig &a) {
return "<FileConfig: " + a.to_string() + ">";
});
py::class_<ScanParameters>(m, "ScanParameters")
.def(py::init<const std::string &>())
.def(py::init<const ScanParameters &>())
.def_property_readonly("enabled", &ScanParameters::enabled)
.def_property_readonly("dac", &ScanParameters::dac)
.def_property_readonly("start", &ScanParameters::start)
.def_property_readonly("stop", &ScanParameters::stop)
.def_property_readonly("step", &ScanParameters::step);
py::class_<ROI>(m, "ROI")
.def(py::init<>())
.def_readwrite("xmin", &ROI::xmin)
.def_readwrite("xmax", &ROI::xmax)
.def_readwrite("ymin", &ROI::ymin)
.def_readwrite("ymax", &ROI::ymax)
.def("__str__", [](const ROI& self){
return fmt::format("ROI: xmin: {} xmax: {} ymin: {} ymax: {}", self.xmin, self.xmax, self.ymin, self.ymax);
})
.def("__repr__", [](const ROI& self){
return fmt::format("<ROI: xmin: {} xmax: {} ymin: {} ymax: {}>", self.xmin, self.xmax, self.ymin, self.ymax);
})
.def("__iter__", [](const ROI &self) {
return py::make_iterator(&self.xmin, &self.ymax+1); //NOLINT
});
py::class_<RawSubFile>(m, "RawSubFile")
.def(py::init<const std::filesystem::path &, DetectorType, size_t,
size_t, size_t>())
.def_property_readonly("bytes_per_frame", &RawSubFile::bytes_per_frame)
.def_property_readonly("pixels_per_frame",
&RawSubFile::pixels_per_frame)
.def("seek", &RawSubFile::seek)
.def("tell", &RawSubFile::tell)
.def_property_readonly("rows", &RawSubFile::rows)
.def_property_readonly("cols", &RawSubFile::cols)
.def("read_frame",
[](RawSubFile &self) {
const uint8_t item_size = self.bytes_per_pixel();
py::array image;
std::vector<ssize_t> shape;
shape.reserve(2);
shape.push_back(self.rows());
shape.push_back(self.cols());
if (item_size == 1) {
image = py::array_t<uint8_t>(shape);
} else if (item_size == 2) {
image = py::array_t<uint16_t>(shape);
} else if (item_size == 4) {
image = py::array_t<uint32_t>(shape);
}
fmt::print("item_size: {} rows: {} cols: {}\n", item_size, self.rows(), self.cols());
self.read_into(
reinterpret_cast<std::byte *>(image.mutable_data()));
return image;
});
// py::class_<ClusterHeader>(m, "ClusterHeader")
// .def(py::init<>())
// .def_readwrite("frame_number", &ClusterHeader::frame_number)
// .def_readwrite("n_clusters", &ClusterHeader::n_clusters)
// .def("__repr__", [](const ClusterHeader &a) { return "<ClusterHeader:
// " + a.to_string() + ">"; });
// py::class_<ClusterV2_>(m, "ClusterV2_")
// .def(py::init<>())
// .def_readwrite("x", &ClusterV2_::x)
// .def_readwrite("y", &ClusterV2_::y)
// .def_readwrite("data", &ClusterV2_::data)
// .def("__repr__", [](const ClusterV2_ &a) { return "<ClusterV2_: " +
// a.to_string(false) + ">"; });
// py::class_<ClusterV2>(m, "ClusterV2")
// .def(py::init<>())
// .def_readwrite("cluster", &ClusterV2::cluster)
// .def_readwrite("frame_number", &ClusterV2::frame_number)
// .def("__repr__", [](const ClusterV2 &a) { return "<ClusterV2: " +
// a.to_string() + ">"; });
// py::class_<ClusterFileV2>(m, "ClusterFileV2")
// .def(py::init<const std::filesystem::path &, const std::string &>())
// .def("read", py::overload_cast<>(&ClusterFileV2::read))
// .def("read", py::overload_cast<int>(&ClusterFileV2::read))
// .def("frame_number", &ClusterFileV2::frame_number)
// .def("write", py::overload_cast<std::vector<ClusterV2> const
// &>(&ClusterFileV2::write))
// .def("close", &ClusterFileV2::close);
// m.def("to_clustV2", [](std::vector<DynamicCluster> &clusters, const int
// frame_number) {
// std::vector<ClusterV2> clusters_;
// for (auto &c : clusters) {
// ClusterV2 cluster;
// cluster.cluster.x = c.x;
// cluster.cluster.y = c.y;
// int i=0;
// for(auto &d : cluster.cluster.data) {
// d=c.get<double>(i++);
// }
// cluster.frame_number = frame_number;
// clusters_.push_back(cluster);
// }
// return clusters_;
// });
}

29
python/src/module.cpp Normal file
View File

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

116
python/src/np_helper.hpp Normal file
View File

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

52
python/src/pedestal.hpp Normal file
View File

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

41
python/src/pixel_map.hpp Normal file
View File

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

106
python/src/raw_file.hpp Normal file
View File

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

View File

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

View File

@ -0,0 +1,46 @@
#include "aare/VarClusterFinder.hpp"
#include "np_helper.hpp"
// #include "aare/defs.hpp"
// #include "aare/fClusterFileV2.hpp"
#include <cstdint>
// #include <filesystem>
#include <pybind11/numpy.h>
// #include <pybind11/iostream.h>
#include <pybind11/pybind11.h>
#include <pybind11/stl.h>
// #include <pybind11/stl/filesystem.h>
// #include <string>
namespace py = pybind11;
using namespace::aare;
void define_var_cluster_finder_bindings(py::module &m) {
PYBIND11_NUMPY_DTYPE(VarClusterFinder<double>::Hit, size, row, col,
reserved, energy, max);
py::class_<VarClusterFinder<double>>(m, "VarClusterFinder")
.def(py::init<Shape<2>, double>())
.def("labeled",
[](VarClusterFinder<double> &self) {
auto ptr = new NDArray<int, 2>(self.labeled());
return return_image_data(ptr);
})
.def("find_clusters",
[](VarClusterFinder<double> &self,
py::array_t<double, py::array::c_style | py::array::forcecast>
img) {
auto view = make_view_2d(img);
self.find_clusters(view);
})
.def("steal_hits",
[](VarClusterFinder<double> &self) {
auto ptr = new std::vector<VarClusterFinder<double>::Hit>(
self.steal_hits());
return return_vector(ptr);
})
.def("total_clusters", &VarClusterFinder<double>::total_clusters);
}

427
src/ClusterFile.cpp Normal file
View File

@ -0,0 +1,427 @@
#include "aare/ClusterFile.hpp"
#include <algorithm>
namespace aare {
ClusterFile::ClusterFile(const std::filesystem::path &fname, size_t chunk_size,
const std::string &mode)
: m_chunk_size(chunk_size), m_mode(mode) {
if (mode == "r") {
fp = fopen(fname.c_str(), "rb");
if (!fp) {
throw std::runtime_error("Could not open file for reading: " +
fname.string());
}
} else if (mode == "w") {
fp = fopen(fname.c_str(), "wb");
if (!fp) {
throw std::runtime_error("Could not open file for writing: " +
fname.string());
}
} else {
throw std::runtime_error("Unsupported mode: " + mode);
}
}
ClusterFile::~ClusterFile() { close(); }
void ClusterFile::close() {
if (fp) {
fclose(fp);
fp = nullptr;
}
}
void ClusterFile::write_frame(int32_t frame_number,
const ClusterVector<int32_t> &clusters) {
if (m_mode != "w") {
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");
}
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.element_offset(), clusters.size(), fp);
// write clusters
// fwrite(clusters.data(), sizeof(Cluster), clusters.size(), fp);
}
std::vector<Cluster3x3> ClusterFile::read_clusters(size_t n_clusters) {
if (m_mode != "r") {
throw std::runtime_error("File not opened for reading");
}
std::vector<Cluster3x3> clusters(n_clusters);
int32_t iframe = 0; // frame number needs to be 4 bytes!
size_t nph_read = 0;
uint32_t nn = m_num_left;
uint32_t nph = m_num_left; // number of clusters in frame needs to be 4
auto buf = reinterpret_cast<Cluster3x3 *>(clusters.data());
// 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(reinterpret_cast<void *>(buf + nph_read),
sizeof(Cluster3x3), 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)) {
// 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(reinterpret_cast<void *>(buf + nph_read),
sizeof(Cluster3x3), 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);
return clusters;
}
std::vector<Cluster3x3> ClusterFile::read_frame(int32_t &out_fnum) {
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");
}
if (fread(&out_fnum, sizeof(out_fnum), 1, fp) != 1) {
throw std::runtime_error("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("Could not read number of clusters");
}
std::vector<Cluster3x3> clusters(n_clusters);
if (fread(clusters.data(), sizeof(Cluster3x3), n_clusters, fp) !=
static_cast<size_t>(n_clusters)) {
throw std::runtime_error("Could not read clusters");
}
return clusters;
}
std::vector<Cluster3x3> ClusterFile::read_cluster_with_cut(size_t n_clusters,
double *noise_map,
int nx, int ny) {
if (m_mode != "r") {
throw std::runtime_error("File not opened for reading");
}
std::vector<Cluster3x3> clusters(n_clusters);
// size_t read_clusters_with_cut(FILE *fp, size_t n_clusters, Cluster *buf,
// uint32_t *n_left, double *noise_map, int
// nx, int ny) {
int iframe = 0;
// uint32_t nph = *n_left;
uint32_t nph = m_num_left;
// uint32_t nn = *n_left;
uint32_t nn = m_num_left;
size_t nph_read = 0;
int32_t t2max, tot1;
int32_t tot3;
// Cluster *ptr = buf;
Cluster3x3 *ptr = clusters.data();
int good = 1;
double noise;
// read photons left from previous frame
if (noise_map)
printf("Using noise map\n");
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;
}
for (size_t iph = 0; iph < nn; iph++) {
// read photons 1 by 1
size_t n_read =
fread(reinterpret_cast<void *>(ptr), sizeof(Cluster3x3), 1, fp);
if (n_read != 1) {
clusters.resize(nph_read);
return clusters;
}
// TODO! error handling on read
good = 1;
if (noise_map) {
if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 && ptr->y < ny) {
tot1 = ptr->data[4];
analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL, NULL, NULL,
NULL);
noise = noise_map[ptr->y * nx + ptr->x];
if (tot1 > noise || t2max > 2 * noise || tot3 > 3 * noise) {
;
} else {
good = 0;
printf("%d %d %f %d %d %d\n", ptr->x, ptr->y, noise,
tot1, t2max, tot3);
}
} else {
printf("Bad pixel number %d %d\n", ptr->x, ptr->y);
good = 0;
}
}
if (good) {
ptr++;
nph_read++;
}
(m_num_left)--;
if (nph_read >= n_clusters)
break;
}
}
if (nph_read < n_clusters) {
// // keep on reading frames and photons until reaching
// n_clusters
while (fread(&iframe, sizeof(iframe), 1, fp)) {
// // printf("%d\n",nph_read);
if (fread(&nph, sizeof(nph), 1, fp)) {
// // printf("** %d\n",nph);
m_num_left = nph;
for (size_t iph = 0; iph < nph; iph++) {
// // read photons 1 by 1
size_t n_read = fread(reinterpret_cast<void *>(ptr),
sizeof(Cluster3x3), 1, fp);
if (n_read != 1) {
clusters.resize(nph_read);
return clusters;
// return nph_read;
}
good = 1;
if (noise_map) {
if (ptr->x >= 0 && ptr->x < nx && ptr->y >= 0 &&
ptr->y < ny) {
tot1 = ptr->data[4];
analyze_cluster(*ptr, &t2max, &tot3, NULL, NULL,
NULL, NULL, NULL);
// noise = noise_map[ptr->y * nx + ptr->x];
noise = noise_map[ptr->y + ny * ptr->x];
if (tot1 > noise || t2max > 2 * noise ||
tot3 > 3 * noise) {
;
} else
good = 0;
} else {
printf("Bad pixel number %d %d\n", ptr->x, ptr->y);
good = 0;
}
}
if (good) {
ptr++;
nph_read++;
}
(m_num_left)--;
if (nph_read >= n_clusters)
break;
}
}
if (nph_read >= n_clusters)
break;
}
}
// printf("%d\n",nph_read);
clusters.resize(nph_read);
return clusters;
}
NDArray<double, 2> calculate_eta2(ClusterVector<int> &clusters) {
NDArray<double, 2> eta2({clusters.size(), 2});
for (size_t i = 0; i < clusters.size(); i++) {
// int32_t t2;
// auto* ptr = reinterpret_cast<int32_t*> (clusters.element_ptr(i) + 2 *
// sizeof(int16_t)); analyze_cluster(clusters.at<Cluster3x3>(i), &t2,
// nullptr, nullptr, &eta2(i,0), &eta2(i,1) , nullptr, nullptr);
auto [x, y] = calculate_eta2(clusters.at<Cluster3x3>(i));
eta2(i, 0) = x;
eta2(i, 1) = y;
}
return eta2;
}
std::array<double, 2> calculate_eta2(Cluster3x3 &cl) {
std::array<double, 2> eta2{};
std::array<int32_t, 4> tot2;
tot2[0] = cl.data[0] + cl.data[1] + cl.data[3] + cl.data[4];
tot2[1] = cl.data[1] + cl.data[2] + cl.data[4] + cl.data[5];
tot2[2] = cl.data[3] + cl.data[4] + cl.data[6] + cl.data[7];
tot2[3] = cl.data[4] + cl.data[5] + cl.data[7] + cl.data[8];
auto c = std::max_element(tot2.begin(), tot2.end()) - tot2.begin();
switch (c) {
case cBottomLeft:
if ((cl.data[3] + cl.data[4]) != 0)
eta2[0] =
static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
if ((cl.data[1] + cl.data[4]) != 0)
eta2[1] =
static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
break;
case cBottomRight:
if ((cl.data[2] + cl.data[5]) != 0)
eta2[0] =
static_cast<double>(cl.data[5]) / (cl.data[4] + cl.data[5]);
if ((cl.data[1] + cl.data[4]) != 0)
eta2[1] =
static_cast<double>(cl.data[4]) / (cl.data[1] + cl.data[4]);
break;
case cTopLeft:
if ((cl.data[7] + cl.data[4]) != 0)
eta2[0] =
static_cast<double>(cl.data[4]) / (cl.data[3] + cl.data[4]);
if ((cl.data[7] + cl.data[4]) != 0)
eta2[1] =
static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
break;
case cTopRight:
if ((cl.data[5] + cl.data[4]) != 0)
eta2[0] =
static_cast<double>(cl.data[5]) / (cl.data[5] + cl.data[4]);
if ((cl.data[7] + cl.data[4]) != 0)
eta2[1] =
static_cast<double>(cl.data[7]) / (cl.data[7] + cl.data[4]);
break;
// default:;
}
return eta2;
}
int analyze_cluster(Cluster3x3 &cl, int32_t *t2, int32_t *t3, char *quad,
double *eta2x, double *eta2y, double *eta3x,
double *eta3y) {
return analyze_data(cl.data, t2, t3, quad, eta2x, eta2y, eta3x, eta3y);
}
int analyze_data(int32_t *data, int32_t *t2, int32_t *t3, char *quad,
double *eta2x, double *eta2y, double *eta3x, double *eta3y) {
int ok = 1;
int32_t tot2[4];
int32_t t2max = 0;
char c = 0;
int32_t val, tot3;
tot3 = 0;
for (int i = 0; i < 4; i++)
tot2[i] = 0;
for (int ix = 0; ix < 3; ix++) {
for (int iy = 0; iy < 3; iy++) {
val = data[iy * 3 + ix];
// printf ("%d ",data[iy * 3 + ix]);
tot3 += val;
if (ix <= 1 && iy <= 1)
tot2[cBottomLeft] += val;
if (ix >= 1 && iy <= 1)
tot2[cBottomRight] += val;
if (ix <= 1 && iy >= 1)
tot2[cTopLeft] += val;
if (ix >= 1 && iy >= 1)
tot2[cTopRight] += val;
}
// printf ("\n");
}
// printf ("\n");
if (t2 || quad) {
t2max = tot2[0];
c = cBottomLeft;
for (int i = 1; i < 4; i++) {
if (tot2[i] > t2max) {
t2max = tot2[i];
c = i;
}
}
// printf("*** %d %d %d %d --
// %d\n",tot2[0],tot2[1],tot2[2],tot2[3],t2max);
if (quad)
*quad = c;
if (t2)
*t2 = t2max;
}
if (t3)
*t3 = tot3;
if (eta2x || eta2y) {
if (eta2x)
*eta2x = 0;
if (eta2y)
*eta2y = 0;
switch (c) {
case cBottomLeft:
if (eta2x && (data[3] + data[4]) != 0)
*eta2x = static_cast<double>(data[4]) / (data[3] + data[4]);
if (eta2y && (data[1] + data[4]) != 0)
*eta2y = static_cast<double>(data[4]) / (data[1] + data[4]);
break;
case cBottomRight:
if (eta2x && (data[2] + data[5]) != 0)
*eta2x = static_cast<double>(data[5]) / (data[4] + data[5]);
if (eta2y && (data[1] + data[4]) != 0)
*eta2y = static_cast<double>(data[4]) / (data[1] + data[4]);
break;
case cTopLeft:
if (eta2x && (data[7] + data[4]) != 0)
*eta2x = static_cast<double>(data[4]) / (data[3] + data[4]);
if (eta2y && (data[7] + data[4]) != 0)
*eta2y = static_cast<double>(data[7]) / (data[7] + data[4]);
break;
case cTopRight:
if (eta2x && t2max != 0)
*eta2x = static_cast<double>(data[5]) / (data[5] + data[4]);
if (eta2y && t2max != 0)
*eta2y = static_cast<double>(data[7]) / (data[7] + data[4]);
break;
default:;
}
}
if (eta3x || eta3y) {
if (eta3x && (data[3] + data[4] + data[5]) != 0)
*eta3x = static_cast<double>(-data[3] + data[3 + 2]) /
(data[3] + data[4] + data[5]);
if (eta3y && (data[1] + data[4] + data[7]) != 0)
*eta3y = static_cast<double>(-data[1] + data[2 * 3 + 1]) /
(data[1] + data[4] + data[7]);
}
return ok;
}
} // namespace aare

View File

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

108
src/ClusterVector.test.cpp Normal file
View File

@ -0,0 +1,108 @@
#include <cstdint>
#include "aare/ClusterVector.hpp"
#include <catch2/matchers/catch_matchers_floating_point.hpp>
#include <catch2/catch_test_macros.hpp>
using aare::ClusterVector;
TEST_CASE("ClusterVector 2x2 int32_t capacity 4, push back then read") {
struct Cluster_i2x2 {
int16_t x;
int16_t y;
int32_t data[4];
};
ClusterVector<int32_t> cv(2, 2, 4);
REQUIRE(cv.capacity() == 4);
REQUIRE(cv.size() == 0);
REQUIRE(cv.cluster_size_x() == 2);
REQUIRE(cv.cluster_size_y() == 2);
// int16_t, int16_t, 2x2 int32_t = 20 bytes
REQUIRE(cv.element_offset() == 20);
//Create a cluster and push back into the vector
Cluster_i2x2 c1 = {1, 2, {3, 4, 5, 6}};
cv.push_back(c1.x, c1.y, reinterpret_cast<std::byte*>(&c1.data[0]));
REQUIRE(cv.size() == 1);
REQUIRE(cv.capacity() == 4);
//Read the cluster back out using copy. TODO! Can we improve the API?
Cluster_i2x2 c2;
std::byte *ptr = cv.element_ptr(0);
std::copy(ptr, ptr + cv.element_offset(), reinterpret_cast<std::byte*>(&c2));
//Check that the data is the same
REQUIRE(c1.x == c2.x);
REQUIRE(c1.y == c2.y);
for(size_t i = 0; i < 4; i++) {
REQUIRE(c1.data[i] == c2.data[i]);
}
}
TEST_CASE("Summing 3x1 clusters of int64"){
struct Cluster_l3x1{
int16_t x;
int16_t y;
int32_t data[3];
};
ClusterVector<int32_t> cv(3, 1, 2);
REQUIRE(cv.capacity() == 2);
REQUIRE(cv.size() == 0);
REQUIRE(cv.cluster_size_x() == 3);
REQUIRE(cv.cluster_size_y() == 1);
//Create a cluster and push back into the vector
Cluster_l3x1 c1 = {1, 2, {3, 4, 5}};
cv.push_back(c1.x, c1.y, reinterpret_cast<std::byte*>(&c1.data[0]));
REQUIRE(cv.capacity() == 2);
REQUIRE(cv.size() == 1);
Cluster_l3x1 c2 = {6, 7, {8, 9, 10}};
cv.push_back(c2.x, c2.y, reinterpret_cast<std::byte*>(&c2.data[0]));
REQUIRE(cv.capacity() == 2);
REQUIRE(cv.size() == 2);
Cluster_l3x1 c3 = {11, 12, {13, 14, 15}};
cv.push_back(c3.x, c3.y, reinterpret_cast<std::byte*>(&c3.data[0]));
REQUIRE(cv.capacity() == 4);
REQUIRE(cv.size() == 3);
auto sums = cv.sum();
REQUIRE(sums.size() == 3);
REQUIRE(sums[0] == 12);
REQUIRE(sums[1] == 27);
REQUIRE(sums[2] == 42);
}
TEST_CASE("Storing floats"){
struct Cluster_f4x2{
int16_t x;
int16_t y;
float data[8];
};
ClusterVector<float> cv(2, 4, 2);
REQUIRE(cv.capacity() == 2);
REQUIRE(cv.size() == 0);
REQUIRE(cv.cluster_size_x() == 2);
REQUIRE(cv.cluster_size_y() == 4);
//Create a cluster and push back into the vector
Cluster_f4x2 c1 = {1, 2, {3.0, 4.0, 5.0, 6.0,3.0, 4.0, 5.0, 6.0}};
cv.push_back(c1.x, c1.y, reinterpret_cast<std::byte*>(&c1.data[0]));
REQUIRE(cv.capacity() == 2);
REQUIRE(cv.size() == 1);
Cluster_f4x2 c2 = {6, 7, {8.0, 9.0, 10.0, 11.0,8.0, 9.0, 10.0, 11.0}};
cv.push_back(c2.x, c2.y, reinterpret_cast<std::byte*>(&c2.data[0]));
REQUIRE(cv.capacity() == 2);
REQUIRE(cv.size() == 2);
auto sums = cv.sum();
REQUIRE(sums.size() == 2);
REQUIRE_THAT(sums[0], Catch::Matchers::WithinAbs(36.0, 1e-6));
REQUIRE_THAT(sums[1], Catch::Matchers::WithinAbs(76.0, 1e-6));
}

74
src/CtbRawFile.cpp Normal file
View File

@ -0,0 +1,74 @@
#include "aare/CtbRawFile.hpp"
#include <fmt/format.h>
namespace aare {
CtbRawFile::CtbRawFile(const std::filesystem::path &fname) : m_master(fname) {
if (m_master.detector_type() != DetectorType::ChipTestBoard) {
throw std::runtime_error(LOCATION + "Not a Ctb file");
}
find_subfiles();
// open the first subfile
m_file.open(m_master.data_fname(0, 0), std::ios::binary);
}
void CtbRawFile::read_into(std::byte *image_buf, DetectorHeader* header) {
if(m_current_frame >= m_master.frames_in_file()){
throw std::runtime_error(LOCATION + " End of file reached");
}
if(m_current_frame != 0 && m_current_frame % m_master.max_frames_per_file() == 0){
open_data_file(m_current_subfile+1);
}
if(header){
m_file.read(reinterpret_cast<char *>(header), sizeof(DetectorHeader));
}else{
m_file.seekg(sizeof(DetectorHeader), std::ios::cur);
}
m_file.read(reinterpret_cast<char *>(image_buf), m_master.image_size_in_bytes());
m_current_frame++;
}
void CtbRawFile::seek(size_t frame_number) {
if (auto index = sub_file_index(frame_number); index != m_current_subfile) {
open_data_file(index);
}
size_t frame_number_in_file = frame_number % m_master.max_frames_per_file();
m_file.seekg((sizeof(DetectorHeader)+m_master.image_size_in_bytes()) * frame_number_in_file);
m_current_frame = frame_number;
}
size_t CtbRawFile::tell() const { return m_current_frame; }
size_t CtbRawFile::image_size_in_bytes() const { return m_master.image_size_in_bytes(); }
size_t CtbRawFile::frames_in_file() const { return m_master.frames_in_file(); }
RawMasterFile CtbRawFile::master() const { return m_master; }
void CtbRawFile::find_subfiles() {
// we can semi safely assume that there is only one module for CTB
while (std::filesystem::exists(m_master.data_fname(0, m_num_subfiles)))
m_num_subfiles++;
fmt::print("Found {} subfiles\n", m_num_subfiles);
}
void CtbRawFile::open_data_file(size_t subfile_index) {
if (subfile_index >= m_num_subfiles) {
throw std::runtime_error(LOCATION + "Subfile index out of range");
}
m_current_subfile = subfile_index;
m_file = std::ifstream(m_master.data_fname(0, subfile_index), std::ios::binary); // only one module for CTB
if (!m_file.is_open()) {
throw std::runtime_error(LOCATION + "Could not open data file");
}
}
} // namespace aare

191
src/Dtype.cpp Normal file
View File

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

54
src/Dtype.test.cpp Normal file
View File

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

79
src/File.cpp Normal file
View File

@ -0,0 +1,79 @@
#include "aare/File.hpp"
#include "aare/NumpyFile.hpp"
#include "aare/RawFile.hpp"
#include <fmt/format.h>
namespace aare {
File::File(const std::filesystem::path &fname, const std::string &mode,
const FileConfig &cfg)
: file_impl(nullptr) {
if (mode != "r") {
throw std::invalid_argument("At the moment only reading is supported");
}
if ((mode == "r") && !std::filesystem::exists(fname)) {
throw std::runtime_error(
fmt::format("File does not exist: {}", fname.string()));
}
// Assuming we are pointing at a master file?
// TODO! How do we read raw files directly?
if (fname.extension() == ".raw" || fname.extension() == ".json") {
// file_impl = new RawFile(fname, mode, cfg);
file_impl = std::make_unique<RawFile>(fname, mode);
}
else if (fname.extension() == ".npy") {
// file_impl = new NumpyFile(fname, mode, cfg);
file_impl = std::make_unique<NumpyFile>(fname, mode, cfg);
} else {
throw std::runtime_error("Unsupported file type");
}
}
File::File(File &&other) noexcept{
std::swap(file_impl, other.file_impl);
}
File& File::operator=(File &&other) noexcept {
if (this != &other) {
File tmp(std::move(other));
std::swap(file_impl, tmp.file_impl);
}
return *this;
}
Frame File::read_frame() { return file_impl->read_frame(); }
Frame File::read_frame(size_t frame_index) {
return file_impl->read_frame(frame_index);
}
size_t File::total_frames() const { return file_impl->total_frames(); }
std::vector<Frame> File::read_n(size_t n_frames) {
return file_impl->read_n(n_frames);
}
void File::read_into(std::byte *image_buf) { file_impl->read_into(image_buf); }
void File::read_into(std::byte *image_buf, size_t n_frames) {
file_impl->read_into(image_buf, n_frames);
}
size_t File::frame_number() { return file_impl->frame_number(tell()); }
size_t File::frame_number(size_t frame_index) {
return file_impl->frame_number(frame_index);
}
size_t File::bytes_per_frame() const { return file_impl->bytes_per_frame(); }
size_t File::pixels_per_frame() const{ return file_impl->pixels_per_frame(); }
void File::seek(size_t frame_index) { file_impl->seek(frame_index); }
size_t File::tell() const { return file_impl->tell(); }
size_t File::rows() const { return file_impl->rows(); }
size_t File::cols() const { return file_impl->cols(); }
size_t File::bitdepth() const { return file_impl->bitdepth(); }
size_t File::bytes_per_pixel() const { return file_impl->bitdepth() / 8; }
DetectorType File::detector_type() const { return file_impl->detector_type(); }
} // namespace aare

74
src/Frame.cpp Normal file
View File

@ -0,0 +1,74 @@
#include "aare/Frame.hpp"
#include <cstddef>
#include <cstring>
#include <iostream>
#include <sys/types.h>
namespace aare {
Frame::Frame(const std::byte *bytes, uint32_t rows, uint32_t cols, Dtype dtype)
: m_rows(rows), m_cols(cols), m_dtype(dtype),
m_data(new std::byte[rows * cols * m_dtype.bytes()]) {
std::memcpy(m_data, bytes, rows * cols * m_dtype.bytes());
}
Frame::Frame(uint32_t rows, uint32_t cols, Dtype dtype)
: m_rows(rows), m_cols(cols), m_dtype(dtype),
m_data(new std::byte[rows * cols * dtype.bytes()]) {
std::memset(m_data, 0, rows * cols * dtype.bytes());
}
uint32_t Frame::rows() const { return m_rows; }
uint32_t Frame::cols() const { return m_cols; }
size_t Frame::bitdepth() const { return m_dtype.bitdepth(); }
Dtype Frame::dtype() const { return m_dtype; }
uint64_t Frame::size() const { return m_rows * m_cols; }
size_t Frame::bytes() const { return m_rows * m_cols * m_dtype.bytes(); }
std::byte *Frame::data() const { return m_data; }
std::byte *Frame::pixel_ptr(uint32_t row, uint32_t col) const{
if ((row >= m_rows) || (col >= m_cols)) {
std::cerr << "Invalid row or column index" << '\n';
return nullptr;
}
return m_data + (row * m_cols + col) * (m_dtype.bytes());
}
Frame &Frame::operator=(Frame &&other) noexcept {
if (this == &other) {
return *this;
}
m_rows = other.rows();
m_cols = other.cols();
m_dtype = other.dtype();
if (m_data != nullptr) {
delete[] m_data;
}
m_data = other.m_data;
other.m_data = nullptr;
other.m_rows = other.m_cols = 0;
other.m_dtype = Dtype(Dtype::TypeIndex::ERROR);
return *this;
}
Frame::Frame(Frame &&other) noexcept
: m_rows(other.rows()), m_cols(other.cols()), m_dtype(other.dtype()),
m_data(other.m_data) {
other.m_data = nullptr;
other.m_rows = other.m_cols = 0;
other.m_dtype = Dtype(Dtype::TypeIndex::ERROR);
}
Frame Frame::clone() const {
Frame frame(m_rows, m_cols, m_dtype);
std::memcpy(frame.m_data, m_data, m_rows * m_cols * m_dtype.bytes());
return frame;
}
} // namespace aare

153
src/Frame.test.cpp Normal file
View File

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

382
src/NDArray.test.cpp Normal file
View File

@ -0,0 +1,382 @@
#include "aare/NDArray.hpp"
#include <array>
#include <catch2/benchmark/catch_benchmark.hpp>
#include <catch2/catch_test_macros.hpp>
using aare::NDArray;
using aare::NDView;
using aare::Shape;
TEST_CASE("Initial size is zero if no size is specified") {
NDArray<double> a;
REQUIRE(a.size() == 0);
REQUIRE(a.shape() == Shape<2>{0, 0});
}
TEST_CASE("Construct from an NDView") {
std::vector<int> some_data(9, 42);
NDView<int, 2> view(some_data.data(), Shape<2>{3, 3});
NDArray<int, 2> image(view);
REQUIRE(image.shape() == view.shape());
REQUIRE(image.size() == view.size());
REQUIRE(image.data() != view.data());
for (uint32_t i = 0; i < image.size(); ++i) {
REQUIRE(image(i) == view(i));
}
// Changing the image doesn't change the view
image = 43;
for (uint32_t i = 0; i < image.size(); ++i) {
REQUIRE(image(i) != view(i));
}
}
TEST_CASE("1D image") {
std::array<int64_t, 1> shape{{20}};
NDArray<short, 1> img(shape, 3);
REQUIRE(img.size() == 20);
REQUIRE(img(5) == 3);
}
TEST_CASE("Accessing a const object") {
const NDArray<double, 3> img({3, 4, 5}, 0);
REQUIRE(img(1, 1, 1) == 0);
REQUIRE(img.size() == 3 * 4 * 5);
REQUIRE(img.shape() == Shape<3>{3, 4, 5});
REQUIRE(img.shape(0) == 3);
REQUIRE(img.shape(1) == 4);
REQUIRE(img.shape(2) == 5);
}
TEST_CASE("Indexing of a 2D image") {
std::array<int64_t, 2> shape{{3, 7}};
NDArray<long> img(shape, 5);
for (uint32_t i = 0; i != img.size(); ++i) {
REQUIRE(img(i) == 5);
}
for (uint32_t i = 0; i != img.size(); ++i) {
img(i) = i;
}
REQUIRE(img(0, 0) == 0);
REQUIRE(img(0, 1) == 1);
REQUIRE(img(1, 0) == 7);
}
TEST_CASE("Indexing of a 3D image") {
NDArray<float, 3> img{{{3, 4, 2}}, 5.0f};
for (uint32_t i = 0; i != img.size(); ++i) {
REQUIRE(img(i) == 5.0f);
}
// Double check general properties
REQUIRE(img.size() == 3 * 4 * 2);
for (uint32_t i = 0; i != img.size(); ++i) {
img(i) = float(i);
}
REQUIRE(img(0, 0, 0) == 0);
REQUIRE(img(0, 0, 1) == 1);
REQUIRE(img(0, 1, 1) == 3);
REQUIRE(img(1, 2, 0) == 12);
REQUIRE(img(2, 3, 1) == 23);
}
TEST_CASE("Divide double by int") {
NDArray<double, 1> a{{5}, 5};
NDArray<int, 1> b{{5}, 5};
a /= b;
for (auto it : a) {
REQUIRE(it == 1.0);
}
}
TEST_CASE("Elementwise multiplication of 3D image") {
std::array<int64_t, 3> shape{3, 4, 2};
NDArray<double, 3> a{shape};
NDArray<double, 3> b{shape};
for (uint32_t i = 0; i != a.size(); ++i) {
a(i) = i;
b(i) = i;
}
// auto c = a * b; // This works but the result is a lazy ArrayMul object
NDArray<double, 3> c = a * b;
REQUIRE(c(0, 0, 0) == 0 * 0);
REQUIRE(c(0, 0, 1) == 1 * 1);
REQUIRE(c(0, 1, 1) == 3 * 3);
REQUIRE(c(1, 2, 0) == 12 * 12);
REQUIRE(c(2, 3, 1) == 23 * 23);
}
NDArray<int> MultiplyNDArrayUsingOperator(NDArray<int> &a, NDArray<int> &b) {
// return a * a * b * b;
NDArray<int>c = a*b;
return c;
}
NDArray<int> MultiplyNDArrayUsingIndex(NDArray<int> &a, NDArray<int> &b) {
NDArray<int> res(a.shape());
for (uint32_t i = 0; i < a.size(); i++) {
// res(i) = a(i) * a(i) * b(i) * b(i);
res(i) = a(i) * b(i);
}
return res;
}
NDArray<int> AddNDArrayUsingOperator(NDArray<int> &a, NDArray<int> &b) {
// return a * a * b * b;
// NDArray<int>c = a+b;
NDArray<int> c(a.shape());
c = a + b;
return c;
}
NDArray<int> AddNDArrayUsingIndex(NDArray<int> &a, NDArray<int> &b) {
NDArray<int> res(a.shape());
for (uint32_t i = 0; i < a.size(); i++) {
// res(i) = a(i) * a(i) * b(i) * b(i);
res(i) = a(i) + b(i);
}
return res;
}
TEST_CASE("Compare two images") {
NDArray<int> a;
NDArray<int> b;
CHECK((a == b));
a = NDArray<int>{{5, 10}, 0};
CHECK((a != b));
b = NDArray<int>{{5, 10}, 0};
CHECK((a == b));
b(3, 3) = 7;
CHECK((a != b));
}
TEST_CASE("Size and shape matches") {
int64_t w = 15;
int64_t h = 75;
std::array<int64_t, 2> shape{w, h};
NDArray<double> a{shape};
REQUIRE(a.size() == static_cast<uint64_t>(w * h));
REQUIRE(a.shape() == shape);
}
TEST_CASE("Initial value matches for all elements") {
double v = 4.35;
NDArray<double> a{{5, 5}, v};
for (uint32_t i = 0; i < a.size(); ++i) {
REQUIRE(a(i) == v);
}
}
TEST_CASE("Data layout of 3D image, fast index last") {
NDArray<int, 3> a{{3, 3, 3}, 0};
REQUIRE(a.size() == 27);
int *ptr = a.data();
for (int i = 0; i < 9; ++i) {
*ptr++ = 10 + i;
REQUIRE(a(0, 0, i) == 10 + i);
REQUIRE(a(i) == 10 + i);
}
}
TEST_CASE("Bitwise and on data") {
NDArray<uint16_t, 1> a({3}, 0);
uint16_t mask = 0x3FF;
a(0) = 16684;
a(1) = 33068;
a(2) = 52608;
a &= mask;
REQUIRE(a(0) == 300);
REQUIRE(a(1) == 300);
REQUIRE(a(2) == 384);
}
TEST_CASE("Elementwise operations on images") {
std::array<int64_t, 2> shape{5, 5};
double a_val = 3.0;
double b_val = 8.0;
SECTION("Add two images") {
NDArray<double> A(shape, a_val);
NDArray<double> B(shape, b_val);
NDArray<double> C = A + B;
// auto C = A+B; // This works but the result is a lazy ArraySum object
// Value of C matches
for (uint32_t i = 0; i < C.size(); ++i) {
REQUIRE(C(i) == a_val + b_val);
}
// Value of A is not changed
for (uint32_t i = 0; i < A.size(); ++i) {
REQUIRE(A(i) == a_val);
}
// Value of B is not changed
for (uint32_t i = 0; i < B.size(); ++i) {
REQUIRE(B(i) == b_val);
}
// A, B and C referes to different data
REQUIRE(A.data() != B.data());
REQUIRE(B.data() != C.data());
}
SECTION("Subtract two images") {
NDArray<double> A(shape, a_val);
NDArray<double> B(shape, b_val);
NDArray<double> C = A - B;
// auto C = A - B; // This works but the result is a lazy ArraySub object
// Value of C matches
for (uint32_t i = 0; i < C.size(); ++i) {
REQUIRE(C(i) == a_val - b_val);
}
// Value of A is not changed
for (uint32_t i = 0; i < A.size(); ++i) {
REQUIRE(A(i) == a_val);
}
// Value of B is not changed
for (uint32_t i = 0; i < B.size(); ++i) {
REQUIRE(B(i) == b_val);
}
// A, B and C referes to different data
REQUIRE(A.data() != B.data());
REQUIRE(B.data() != C.data());
}
SECTION("Multiply two images") {
NDArray<double> A(shape, a_val);
NDArray<double> B(shape, b_val);
// auto C = A * B; // This works but the result is a lazy ArrayMul object
NDArray<double> C = A * B;
// Value of C matches
for (uint32_t i = 0; i < C.size(); ++i) {
REQUIRE(C(i) == a_val * b_val);
}
// Value of A is not changed
for (uint32_t i = 0; i < A.size(); ++i) {
REQUIRE(A(i) == a_val);
}
// Value of B is not changed
for (uint32_t i = 0; i < B.size(); ++i) {
REQUIRE(B(i) == b_val);
}
// A, B and C referes to different data
REQUIRE(A.data() != B.data());
REQUIRE(B.data() != C.data());
}
SECTION("Divide two images") {
NDArray<double> A(shape, a_val);
NDArray<double> B(shape, b_val);
// auto C = A / B; // This works but the result is a lazy ArrayDiv object
NDArray<double> C = A / B;
// Value of C matches
for (uint32_t i = 0; i < C.size(); ++i) {
REQUIRE(C(i) == a_val / b_val);
}
// Value of A is not changed
for (uint32_t i = 0; i < A.size(); ++i) {
REQUIRE(A(i) == a_val);
}
// Value of B is not changed
for (uint32_t i = 0; i < B.size(); ++i) {
REQUIRE(B(i) == b_val);
}
// A, B and C referes to different data
REQUIRE(A.data() != B.data());
REQUIRE(B.data() != C.data());
}
SECTION("subtract scalar") {
NDArray<double> A(shape, a_val);
NDArray<double> B(shape, b_val);
double v = 1.0;
auto C = A - v;
REQUIRE(C.data() != A.data());
// Value of C matches
for (uint32_t i = 0; i < C.size(); ++i) {
REQUIRE(C(i) == a_val - v);
}
// Value of A is not changed
for (uint32_t i = 0; i < A.size(); ++i) {
REQUIRE(A(i) == a_val);
}
}
SECTION("add scalar") {
NDArray<double> A(shape, a_val);
NDArray<double> B(shape, b_val);
double v = 1.0;
auto C = A + v;
REQUIRE(C.data() != A.data());
// Value of C matches
for (uint32_t i = 0; i < C.size(); ++i) {
REQUIRE(C(i) == a_val + v);
}
// Value of A is not changed
for (uint32_t i = 0; i < A.size(); ++i) {
REQUIRE(A(i) == a_val);
}
}
SECTION("divide with scalar") {
NDArray<double> A(shape, a_val);
NDArray<double> B(shape, b_val);
double v = 3.7;
auto C = A / v;
REQUIRE(C.data() != A.data());
// Value of C matches
for (uint32_t i = 0; i < C.size(); ++i) {
REQUIRE(C(i) == a_val / v);
}
// Value of A is not changed
for (uint32_t i = 0; i < A.size(); ++i) {
REQUIRE(A(i) == a_val);
}
}
SECTION("multiply with scalar") {
NDArray<double> A(shape, a_val);
NDArray<double> B(shape, b_val);
double v = 3.7;
auto C = A / v;
REQUIRE(C.data() != A.data());
// Value of C matches
for (uint32_t i = 0; i < C.size(); ++i) {
REQUIRE(C(i) == a_val / v);
}
// Value of A is not changed
for (uint32_t i = 0; i < A.size(); ++i) {
REQUIRE(A(i) == a_val);
}
}
}

193
src/NDView.test.cpp Normal file
View File

@ -0,0 +1,193 @@
#include "aare/NDView.hpp"
#include <catch2/catch_test_macros.hpp>
#include <iostream>
#include <vector>
using aare::NDView;
using aare::Shape;
TEST_CASE("Element reference 1D") {
std::vector<int> vec;
for (int i = 0; i != 10; ++i) {
vec.push_back(i);
}
NDView<int, 1> data(vec.data(), Shape<1>{10});
REQUIRE(vec.size() == static_cast<size_t>(data.size()));
for (int i = 0; i != 10; ++i) {
REQUIRE(data(i) == vec[i]);
REQUIRE(data[i] == vec[i]);
}
}
TEST_CASE("Element reference 2D") {
std::vector<int> vec;
for (int i = 0; i != 12; ++i) {
vec.push_back(i);
}
NDView<int, 2> data(vec.data(), Shape<2>{3, 4});
REQUIRE(vec.size() == static_cast<size_t>(data.size()));
int i = 0;
for (int row = 0; row != 3; ++row) {
for (int col = 0; col != 4; ++col) {
REQUIRE(data(row, col) == i);
REQUIRE(data[i] == vec[i]);
++i;
}
}
}
TEST_CASE("Element reference 3D") {
std::vector<int> vec;
for (int i = 0; i != 24; ++i) {
vec.push_back(i);
}
NDView<int, 3> data(vec.data(), Shape<3>{2, 3, 4});
REQUIRE(vec.size() == static_cast<size_t>(data.size()));
int i = 0;
for (int frame = 0; frame != 2; ++frame) {
for (int row = 0; row != 3; ++row) {
for (int col = 0; col != 4; ++col) {
REQUIRE(data(frame, row, col) == i);
REQUIRE(data[i] == vec[i]);
++i;
}
}
}
}
TEST_CASE("Plus and miuns with single value") {
std::vector<int> vec;
for (int i = 0; i != 12; ++i) {
vec.push_back(i);
}
NDView<int, 2> data(vec.data(), Shape<2>{3, 4});
data += 5;
int i = 0;
for (int row = 0; row != 3; ++row) {
for (int col = 0; col != 4; ++col) {
REQUIRE(data(row, col) == i + 5);
++i;
}
}
data -= 3;
i = 0;
for (int row = 0; row != 3; ++row) {
for (int col = 0; col != 4; ++col) {
REQUIRE(data(row, col) == i + 2);
++i;
}
}
}
TEST_CASE("Multiply and divide with single value") {
std::vector<int> vec;
for (int i = 0; i != 12; ++i) {
vec.push_back(i);
}
NDView<int, 2> data(vec.data(), Shape<2>{3, 4});
data *= 5;
int i = 0;
for (int row = 0; row != 3; ++row) {
for (int col = 0; col != 4; ++col) {
REQUIRE(data(row, col) == i * 5);
++i;
}
}
data /= 3;
i = 0;
for (int row = 0; row != 3; ++row) {
for (int col = 0; col != 4; ++col) {
REQUIRE(data(row, col) == (i * 5) / 3);
++i;
}
}
}
TEST_CASE("elementwise assign") {
std::vector<int> vec(25);
NDView<int, 2> data(vec.data(), Shape<2>{5, 5});
data = 3;
for (auto it : data) {
REQUIRE(it == 3);
}
}
TEST_CASE("iterators") {
std::vector<int> vec;
for (int i = 0; i != 12; ++i) {
vec.push_back(i);
}
NDView<int, 1> data(vec.data(), Shape<1>{12});
int i = 0;
for (const auto item : data) {
REQUIRE(item == vec[i]);
++i;
}
REQUIRE(i == 12);
for (auto ptr = data.begin(); ptr != data.end(); ++ptr) {
*ptr += 1;
}
for (auto &item : data) {
++item;
}
i = 0;
for (const auto item : data) {
REQUIRE(item == i + 2);
++i;
}
}
// TEST_CASE("shape from vector") {
// std::vector<int> vec;
// for (int i = 0; i != 12; ++i) {
// vec.push_back(i);
// }
// std::vector<int64_t> shape{3, 4};
// NDView<int, 2> data(vec.data(), shape);
// }
TEST_CASE("divide with another span") {
std::vector<int> vec0{9, 12, 3};
std::vector<int> vec1{3, 2, 1};
std::vector<int> result{3, 6, 3};
NDView<int, 1> data0(vec0.data(), Shape<1>{static_cast<int64_t>(vec0.size())});
NDView<int, 1> data1(vec1.data(), Shape<1>{static_cast<int64_t>(vec1.size())});
data0 /= data1;
for (size_t i = 0; i != vec0.size(); ++i) {
REQUIRE(data0[i] == result[i]);
}
}
TEST_CASE("Retrieve shape") {
std::vector<int> vec;
for (int i = 0; i != 12; ++i) {
vec.push_back(i);
}
NDView<int, 2> data(vec.data(), Shape<2>{3, 4});
REQUIRE(data.shape()[0] == 3);
REQUIRE(data.shape()[1] == 4);
}
TEST_CASE("compare two views") {
std::vector<int> vec1;
for (int i = 0; i != 12; ++i) {
vec1.push_back(i);
}
NDView<int, 2> view1(vec1.data(), Shape<2>{3, 4});
std::vector<int> vec2;
for (int i = 0; i != 12; ++i) {
vec2.push_back(i);
}
NDView<int, 2> view2(vec2.data(), Shape<2>{3, 4});
REQUIRE((view1 == view2));
}

200
src/NumpyFile.cpp Normal file
View File

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

50
src/NumpyFile.test.cpp Normal file
View File

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

269
src/NumpyHelpers.cpp Normal file
View File

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

62
src/NumpyHelpers.test.cpp Normal file
View File

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

103
src/Pedestal.test.cpp Normal file
View File

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

130
src/PixelMap.cpp Normal file
View File

@ -0,0 +1,130 @@
#include "aare/PixelMap.hpp"
#include <array>
namespace aare {
NDArray<ssize_t, 2> GenerateMoench03PixelMap() {
std::array<int, 32> const adc_nr = {300, 325, 350, 375, 300, 325, 350, 375,
200, 225, 250, 275, 200, 225, 250, 275,
100, 125, 150, 175, 100, 125, 150, 175,
0, 25, 50, 75, 0, 25, 50, 75};
int const sc_width = 25;
int const nadc = 32;
int const pixels_per_sc = 5000;
NDArray<ssize_t, 2> order_map({400, 400});
int pixel = 0;
for (int i = 0; i != pixels_per_sc; ++i) {
for (int i_adc = 0; i_adc != nadc; ++i_adc) {
int const col = adc_nr[i_adc] + (i % sc_width);
int row = 0;
if ((i_adc / 4) % 2 == 0)
row = 199 - (i / sc_width);
else
row = 200 + (i / sc_width);
order_map(row, col) = pixel;
pixel++;
}
}
return order_map;
}
NDArray<ssize_t, 2> GenerateMoench05PixelMap() {
std::array<int, 3> adc_numbers = {5, 9, 1};
NDArray<ssize_t, 2> order_map({160, 150});
int n_pixel = 0;
for (int row = 0; row < 160; row++) {
for (int i_col = 0; i_col < 50; i_col++) {
n_pixel = row * 50 + i_col;
for (int i_sc = 0; i_sc < 3; i_sc++) {
int col = 50 * i_sc + i_col;
int adc_nr = adc_numbers[i_sc];
int i_analog = n_pixel * 12 + adc_nr;
// analog_frame[row * 150 + col] = analog_data[i_analog] & 0x3FFF;
order_map(row, col) = i_analog;
}
}
}
return order_map;
}
NDArray<ssize_t, 2> GenerateMoench05PixelMap1g() {
std::array<int, 3> adc_numbers = {1, 2, 0};
NDArray<ssize_t, 2> order_map({160, 150});
int n_pixel = 0;
for (int row = 0; row < 160; row++) {
for (int i_col = 0; i_col < 50; i_col++) {
n_pixel = row * 50 + i_col;
for (int i_sc = 0; i_sc < 3; i_sc++) {
int col = 50 * i_sc + i_col;
int adc_nr = adc_numbers[i_sc];
int i_analog = n_pixel * 3 + adc_nr;
// analog_frame[row * 150 + col] = analog_data[i_analog] & 0x3FFF;
order_map(row, col) = i_analog;
}
}
}
return order_map;
}
NDArray<ssize_t, 2> GenerateMoench05PixelMapOld() {
std::array<int, 3> adc_numbers = {9, 13, 1};
NDArray<ssize_t, 2> order_map({160, 150});
int n_pixel = 0;
for (int row = 0; row < 160; row++) {
for (int i_col = 0; i_col < 50; i_col++) {
n_pixel = row * 50 + i_col;
for (int i_sc = 0; i_sc < 3; i_sc++) {
int col = 50 * i_sc + i_col;
int adc_nr = adc_numbers[i_sc];
int i_analog = n_pixel * 32 + adc_nr;
// analog_frame[row * 150 + col] = analog_data[i_analog] & 0x3FFF;
order_map(row, col) = i_analog;
}
}
}
return order_map;
}
NDArray<ssize_t, 2>GenerateEigerFlipRowsPixelMap(){
NDArray<ssize_t, 2> order_map({256, 512});
for(int row = 0; row < 256; row++){
for(int col = 0; col < 512; col++){
order_map(row, col) = 255*512-row*512 + col;
}
}
return order_map;
}
NDArray<ssize_t, 2>GenerateMH02SingleCounterPixelMap(){
NDArray<ssize_t, 2> order_map({48, 48});
for(int row = 0; row < 48; row++){
for(int col = 0; col < 48; col++){
order_map(row, col) = row*48 + col;
}
}
return order_map;
}
NDArray<ssize_t, 3> GenerateMH02FourCounterPixelMap(){
NDArray<ssize_t, 3> order_map({4, 48, 48});
for (int counter=0; counter<4; counter++){
for(int row = 0; row < 48; row++){
for(int col = 0; col < 48; col++){
order_map(counter, row, col) = counter*48*48 + row*48 + col;
}
}
}
return order_map;
}
} // namespace aare

Some files were not shown because too many files have changed in this diff Show More