diff --git a/.github/workflows/build-test.yml b/.github/workflows/build-test.yml index 25a17f42e4..bc3ae25b03 100644 --- a/.github/workflows/build-test.yml +++ b/.github/workflows/build-test.yml @@ -33,6 +33,9 @@ on: default: false jobs: build: + env: + YARDL_VERSION: 0.6.7 + defaults: {run: {shell: 'bash -leo pipefail {0}'}} runs-on: ${{ matrix.os }} strategy: @@ -44,9 +47,7 @@ jobs: cuda_version: "0" BUILD_FLAGS: "-DSTIR_OPENMP=ON" BUILD_TYPE: "Release" - parallelproj: "ON" - ROOT: "ON" - ITK: "OFF" + PETSIRD: "ON" - os: ubuntu-24.04 compiler: gcc # compiler_version: 9 @@ -56,6 +57,7 @@ jobs: parallelproj: "ON" ROOT: "OFF" ITK: "OFF" + PETSIRD: "OFF" - os: ubuntu-24.04 compiler: clang #compiler_version: @@ -66,6 +68,7 @@ jobs: ROOT: "OFF" # currently using ITK 5.2 which doesn't like clang 14 ITK: "OFF" + PETSIRD: "OFF" - os: ubuntu-24.04 compiler: gcc compiler_version: 10 @@ -75,6 +78,7 @@ jobs: parallelproj: "OFF" ROOT: "OFF" ITK: "ON" + PETSIRD: "OFF" # gcc-12, C++20 test (Interesting, as gcc-12 does not support all of C++20 yet). # However, commented out, as this seems to hang during install step for unknown reasons. # See https://github.com/UCL/STIR/pull/1605 @@ -97,6 +101,7 @@ jobs: ROOT: "OFF" # Currently disabled due to out of disk space, see https://github.com/UCL/STIR/issues/1618 ITK: "OFF" + PETSIRD: "OFF" - os: ubuntu-24.04 # shared library build compiler: gcc @@ -107,6 +112,7 @@ jobs: parallelproj: "ON" ROOT: "OFF" ITK: "ON" + PETSIRD: "OFF" - os: ubuntu-24.04 compiler: gcc # currently CUDA doesn't support gcc 14 yet @@ -117,6 +123,7 @@ jobs: parallelproj: "ON" ROOT: "OFF" ITK: "ON" + PETSIRD: "OFF" - os: macOS-latest compiler: gcc # compiler_version: 11 @@ -127,15 +134,16 @@ jobs: BUILD_TYPE: "Debug" ROOT: "OFF" ITK: "OFF" + PETSIRD: "OFF" - os: macOS-latest - compiler: clang - compiler_version: 18 + compiler: conda-compiler + compiler_version: "" cuda_version: "0" - BUILD_FLAGS: "-DSTIR_OPENMP=OFF" - parallelproj: "OFF" + parallelproj: "ON" BUILD_TYPE: "Release" - ROOT: "OFF" - ITK: "OFF" + ROOT: "ON" + ITK: "ON" + PETSIRD: "ON" - os: macOS-latest compiler: clang compiler_version: 21 @@ -147,35 +155,16 @@ jobs: BUILD_TYPE: "Release" ROOT: "OFF" ITK: "OFF" + PETSIRD: "OFF" # let's run all of them, as opposed to aborting when one fails - fail-fast: false + # TODO remove comment + # fail-fast: false - name: ${{ matrix.os }}-${{ matrix.compiler }}${{ matrix.compiler_version }}-cuda${{ matrix.cuda_version }}-${{ matrix.BUILD_TYPE }}-pp=${{ matrix.parallelproj }}-ROOT=${{ matrix.ROOT }} + name: ${{ matrix.os }}-${{ matrix.compiler }}${{ matrix.compiler_version }}-cuda${{ matrix.cuda_version }}-${{ matrix.BUILD_TYPE }}-petsird=${{ matrix.PETSIRD}}-pp=${{ matrix.parallelproj }}-ROOT=${{ matrix.ROOT }} steps: - - name: disk space - shell: bash - run: | - case ${{matrix.os}} in - (ubuntu* | macOS*) - df -h - # saves about 2GB - echo removing dotnet - sudo rm -rf /usr/share/dotnet - # saves about 10 GB - echo removing agent_tools - sudo rm -rf "$AGENT_TOOLSDIRECTORY" - # saves about 10 GB - echo removing android files - sudo rm -rf /usr/local/lib/android - df -h - ;; - (windows*) - # no idea how to do this - ;; - esac - - uses: actions/checkout@v4 + - uses: actions/checkout@v7 with: submodules: recursive @@ -184,15 +173,93 @@ jobs: run: | case ${{matrix.os}} in (ubuntu* | macOS*) - sudo .github/workflows/GHA_increase_disk_space.sh + if [ "${ACT:-false}" = "false" ]; then + sudo .github/workflows/GHA_increase_disk_space.sh + fi ;; (windows*) # no idea what to do here ;; esac - - name: set_compiler_variables - shell: bash + - if: matrix.cuda_version != '0' + uses: Jimver/cuda-toolkit@v0.2.35 + id: cuda-toolkit + with: + cuda: ${{ matrix.cuda_version }} + linux-local-args: '["--toolkit"]' + + - name: ccache + uses: hendrikmuhs/ccache-action@v1 + with: + key: ${{ matrix.os }}-${{ matrix.compiler }}${{ matrix.compiler_version }}-cuda${{ matrix.cuda_version }}-${{ matrix.BUILD_TYPE }} + max-size: "2G" + + - name: Setup conda environment + if: matrix.PETSIRD == 'ON' + uses: mamba-org/setup-micromamba@v3 + with: + environment-file: build_environment.yml + environment-name: stir + cache-environment: true + cache-downloads: true + + - name: General env variables + run: | + if [[ '${{ matrix.os }}' == macOS* ]]; then + if test ! -z "$CONDA_PREFIX"; then + # see https://conda-forge.org/docs/maintainer/knowledge_base/#newer-c-features-with-old-sdk + echo MACOS: $CXXFLAGS + echo CXXFLAGS=\"${CXXFLAGS} -D_LIBCPP_DISABLE_AVAILABILITY\">> $GITHUB_ENV + echo MACOS: $CXXFLAGS + fi + fi + # We will install all software, including some external dependencies, here + echo CMAKE_INSTALL_PREFIX="${GITHUB_WORKSPACE}/install" >> $GITHUB_ENV + + - name: Checkout PETSIRD + if: matrix.PETSIRD == 'ON' + uses: actions/checkout@v7 + with: + repository: ETSInitiative/PETSIRD + ref: v0.9.1 + path: external/PETSIRD + + - name: Install yardl + if: matrix.PETSIRD == 'ON' + run: | + OS="$(uname -s)" + ARCH="$(uname -m)" + echo $OS + echo $ARCH + case "$OS" in + Linux) + YARDL_OS="linux" + ;; + Darwin) + YARDL_OS="darwin" + ;; + Windows) + YARDL_OS="windows" + ;; + esac + + curl -L https://github.com/microsoft/yardl/releases/download/v${YARDL_VERSION}/yardl_${YARDL_VERSION}_${YARDL_OS}_${ARCH}.tar.gz -o yardl_${YARDL_VERSION}_${YARDL_OS}_${ARCH}.tar.gz + mkdir yardl + tar -xzf "yardl_${YARDL_VERSION}_${YARDL_OS}_${ARCH}.tar.gz" -C yardl + rm "yardl_${YARDL_VERSION}_${YARDL_OS}_${ARCH}.tar.gz" + echo "$PWD/yardl" >> "$GITHUB_PATH" + + - name: Install PETSIRD + if: matrix.PETSIRD == 'ON' + run: | + set -e + micromamba install conda-forge::just + cd external/PETSIRD/ + just cmake_install_prefix="${CMAKE_INSTALL_PREFIX}" build-cpp + + - name: set_compiler_variables (no conda) + if: matrix.PETSIRD !='ON' run: | set -ex if test 'XX${{ matrix.compiler }}' = 'XXclang'; then @@ -229,33 +296,26 @@ jobs: echo CC="$CC" >> $GITHUB_ENV echo CXX="$CXX" >> $GITHUB_ENV - - if: matrix.cuda_version != '0' - uses: Jimver/cuda-toolkit@v0.2.19 - id: cuda-toolkit - with: - cuda: ${{ matrix.cuda_version }} - linux-local-args: '["--toolkit"]' - - - name: install_dependencies - shell: bash + - name: install_dependencies (no conda) + if: matrix.PETSIRD !='ON' run: | set -ex - # We will install some external dependencies here - CMAKE_INSTALL_PREFIX=${GITHUB_WORKSPACE}/install + case ${{matrix.os}} in (ubuntu*) sudo apt update # install compiler if test 'XX${{ matrix.compiler }}' = 'XXclang'; then # package is called clang, need libomp-dev for OpenMP support - sudo apt install $CC libomp-dev + sudo apt install -y $CC libomp-dev else - sudo apt install $CXX + sudo apt install -y $CXX fi + # other dependencies - sudo apt install libboost-dev libhdf5-serial-dev swig python3-dev nlohmann-json3-dev ninja-build + sudo apt install -y cmake libboost-dev libhdf5-serial-dev swig python3-dev nlohmann-json3-dev ninja-build if test "${{matrix.ITK}}XX" == "ONXX"; then - sudo apt install libinsighttoolkit5-dev + sudo apt install -y libinsighttoolkit5-dev fi # free up some disk space apt autoremove --purge && sudo apt clean @@ -268,12 +328,11 @@ jobs: if ! command -v swig > /dev/null; then brew install swig fi - brew install boost nlohmann-json # Temp fix to 3.13 due to https://github.com/UCL/STIR/issues/1638 #if ! command -v python3 > /dev/null; then # brew install python #fi - brew install python@3.13 + brew install boost nlohmann-json python@3.13 brew ls python@3.13 PYTHON_EXECUTABLE=$(find /opt/homebrew/Cellar/python@3.13 -type f -name python3.13) #PYTHON_EXECUTABLE=$(which python3) @@ -286,10 +345,11 @@ jobs: PYTHON_EXECUTABLE=$(which python3) ;; esac - echo PYTHON_EXECUTABLE="$PYTHON_EXECUTABLE" >> $GITHUB_ENV + ${PYTHON_EXECUTABLE} -m venv ${GITHUB_WORKSPACE}/my-env source ${GITHUB_WORKSPACE}/my-env/bin/activate - #python -m pip install -U pip + python -m pip install -U pip + case ${{matrix.os}} in (macOS*) # attempt to get round buggy Accelerate builds, see https://github.com/numpy/numpy/issues/15947 @@ -325,7 +385,7 @@ jobs: if test "${{matrix.ROOT}}XX" == "ONXX"; then case ${{matrix.os}} in (ubuntu*) - sudo apt install libtbb-dev libvdt-dev libgif-dev + sudo apt install -y libtbb-dev libvdt-dev libgif-dev libxpm-dev ROOT_file=root_v6.34.00.Linux-ubuntu24.04-x86_64-gcc13.2.tar.gz #root_v6.34.00.Linux-ubuntu24.10-x86_64-gcc14.2.tar.gz ;; @@ -337,18 +397,20 @@ jobs: tar -xzvf "$ROOT_file" rm "$ROOT_file" source root/bin/thisroot.sh - echo ROOTSYS="$ROOTSYS" >> $GITHUB_ENV + # thisroot.sh adds ROOTSYS to CMAKE_PREFIX_PATH. + echo ROOTSYS="$ROOTSYS" >> "$GITHUB_ENV" + echo CMAKE_PREFIX_PATH="$CMAKE_PREFIX_PATH" >> "$GITHUB_ENV" fi # we'll install some dependencies with shared libraries, so need to let the OS know # thisroot.sh also modified the path, so save that for the recon_test_pack case ${{matrix.os}} in (ubuntu*) - echo LD_LIBRARY_PATH="${CMAKE_INSTALL_PREFIX}/lib:$LD_LIBRARY_PATH" >> $GITHUB_ENV + echo LD_LIBRARY_PATH="${CMAKE_INSTALL_PREFIX}/lib:${LD_LIBRARY_PATH}" >> "$GITHUB_ENV" echo PATH="$PATH" >> $GITHUB_ENV ;; (macOS*) - echo DYLD_FALLBACK_LIBRARY_PATH="${CMAKE_INSTALL_PREFIX}/lib:$DYLD_FALLBACK_LIBRARY_PATH" >> $GITHUB_ENV + echo DYLD_FALLBACK_LIBRARY_PATH="${CMAKE_INSTALL_PREFIX}/lib:${DYLD_FALLBACK_LIBRARY_PATH}" >> "$GITHUB_ENV" echo PATH="$PATH" >> $GITHUB_ENV ;; (windows*) @@ -356,36 +418,35 @@ jobs: ;; esac - - name: ccache - uses: hendrikmuhs/ccache-action@v1 - with: - key: ${{ matrix.os }}-${{ matrix.compiler }}${{ matrix.compiler_version }}-cuda${{ matrix.cuda_version }}-${{ matrix.BUILD_TYPE }} - max-size: "2G" - - name: configure - shell: bash env: BUILD_TYPE: ${{ matrix.BUILD_TYPE }} BUILD_FLAGS: ${{ matrix.BUILD_FLAGS }} run: | set -ex - source ${GITHUB_WORKSPACE}/my-env/bin/activate + if test -f "${GITHUB_WORKSPACE}/my-env/bin/activate"; then + source ${GITHUB_WORKSPACE}/my-env/bin/activate + fi #export PATH="/usr/lib/ccache:/usr/local/opt/ccache/libexec:$PATH" cmake --version if test "XX$CC" != "XX"; then $CC --version $CXX --version fi - CMAKE_INSTALL_PREFIX=${GITHUB_WORKSPACE}/install - # make available to jobs below - echo CMAKE_INSTALL_PREFIX="$CMAKE_INSTALL_PREFIX" >> $GITHUB_ENV if [ -n "$ROOTSYS" ]; then - # make sure we find ROOT (and vdt, which is installed in the same place) - EXTRA_BUILD_FLAGS=-DCMAKE_PREFIX_PATH:PATH="$ROOTSYS" + if test "${{matrix.PETSIRD}}XX" == "ONXX"; then + EXTRA_BUILD_FLAGS="${EXTRA_BUILD_FLAGS} -DCMAKE_PREFIX_PATH=${CONDA_PREFIX}" + EXTRA_BUILD_FLAGS="${EXTRA_BUILD_FLAGS} -DHDF5_ROOT=${CONDA_PREFIX}" + #EXTRA_BUILD_FLAGS="${EXTRA_BUILD_FLAGS} -DVDT_LIBRARY=${CONDA_PREFIX}/lib/libvdt.so" + #EXTRA_BUILD_FLAGS="${EXTRA_BUILD_FLAGS} -DVDT_INCLUDE_DIR=${CONDA_PREFIX}/include" + else + EXTRA_BUILD_FLAGS="${EXTRA_BUILD_FLAGS} -DCMAKE_PREFIX_PATH=${ROOTSYS};/usr" + fi fi EXTRA_BUILD_FLAGS="${EXTRA_BUILD_FLAGS} -DBUILD_SWIG_PYTHON=ON -DPython_EXECUTABLE=`which python`" EXTRA_BUILD_FLAGS="${EXTRA_BUILD_FLAGS} -DCMAKE_INSTALL_PREFIX=${CMAKE_INSTALL_PREFIX} -DCMAKE_BUILD_TYPE=${BUILD_TYPE}" EXTRA_BUILD_FLAGS="${EXTRA_BUILD_FLAGS} -DDOWNLOAD_ZENODO_TEST_DATA=ON" + EXTRA_BUILD_FLAGS="${EXTRA_BUILD_FLAGS} -DGRAPHICS=None" EXTRA_BUILD_FLAGS="${EXTRA_BUILD_FLAGS} -DDISABLE_STIR_LOCAL=OFF -DSTIR_LOCAL=${GITHUB_WORKSPACE}/examples/C++/using_STIR_LOCAL" echo "cmake flags $BUILD_FLAGS $EXTRA_BUILD_FLAGS" mkdir build @@ -409,16 +470,16 @@ jobs: if: ${{ github.event_name == 'workflow_dispatch' && inputs.debug_enabled == 'true' }} - name: build - shell: bash env: BUILD_TYPE: ${{ matrix.BUILD_TYPE }} run: | cd ${GITHUB_WORKSPACE}/build; - source ${GITHUB_WORKSPACE}/my-env/bin/activate - cmake --build . --config ${BUILD_TYPE}} --target install + if test -f "${GITHUB_WORKSPACE}/my-env/bin/activate"; then + source "${GITHUB_WORKSPACE}/my-env/bin/activate" + fi + cmake --build . --config ${BUILD_TYPE} --target install - name: ctest - shell: bash env: BUILD_TYPE: ${{ matrix.BUILD_TYPE }} run: | @@ -448,7 +509,7 @@ jobs: uses: actions/upload-artifact@v4 if: failure() with: - name: ctest_log_files-${{ matrix.os }}-${{ matrix.compiler }}${{ matrix.compiler_version }}-cuda${{ matrix.cuda_version }}-${{ matrix.BUILD_TYPE }}-pp=${{ matrix.parallelproj }}-ROOT=${{ matrix.ROOT }} + name: ctest_log_files-${{ matrix.os }}-${{ matrix.compiler }}${{ matrix.compiler_version }}-cuda${{ matrix.cuda_version }}-${{ matrix.BUILD_TYPE }}-petsird=${{ matrix.PETSIRD}}-pp=${{ matrix.parallelproj }}-ROOT=${{ matrix.ROOT }} path: | ${{ github.workspace }}/build/**/*.log ${{ github.workspace }}/build/**/*.hv @@ -458,7 +519,6 @@ jobs: retention-days: 7 - name: C++ examples with STIR_LOCAL - shell: bash run: | set -ex; PATH=${CMAKE_INSTALL_PREFIX}/bin:$PATH @@ -479,9 +539,7 @@ jobs: ${EXE_LOC}/demo4_obj_fun demo4_obj_fun.par ${EXE_LOC}/demo5_line_search demo5_line_search.par - - name: remove build - shell: bash - # remove to create some disk space + - name: remove build # remove to create some disk space run: | set -x cd ${GITHUB_WORKSPACE} @@ -532,7 +590,7 @@ jobs: uses: actions/upload-artifact@v4 if: failure() with: - name: recon_test_pack_log_files-${{ matrix.os }}-${{ matrix.compiler }}${{ matrix.compiler_version }}-cuda${{ matrix.cuda_version }}-${{ matrix.BUILD_TYPE }}-pp=${{ matrix.parallelproj }}-ROOT=${{ matrix.ROOT }} + name: recon_test_pack_log_files-${{ matrix.os }}-${{ matrix.compiler }}${{ matrix.compiler_version }}-cuda${{ matrix.cuda_version }}-${{ matrix.BUILD_TYPE }}-petsird=${{ matrix.PETSIRD}}-pp=${{ matrix.parallelproj }}-ROOT=${{ matrix.ROOT }} path: | ${{ github.workspace }}/recon_test_pack/**/*.log ${{ github.workspace }}/recon_test_pack/**/my_*v @@ -553,7 +611,9 @@ jobs: shell: bash run: | set -ex - source ${GITHUB_WORKSPACE}/my-env/bin/activate + if test -f "${GITHUB_WORKSPACE}/my-env/bin/activate"; then + source "${GITHUB_WORKSPACE}/my-env/bin/activate" + fi # Run Python tests, making sure we're using the correct Python interpreter which python export PYTHONPATH=${CMAKE_INSTALL_PREFIX}/python diff --git a/CMakeLists.txt b/CMakeLists.txt index 3ae4956c1f..66943261fa 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -109,6 +109,7 @@ option(STIR_ENABLE_EXPERIMENTAL "disable use of STIR experimental code" OFF) # d option(DISABLE_NiftyPET_PROJECTOR "disable use of NiftyPET projector" OFF) option(DISABLE_Parallelproj_PROJECTOR "disable use of Parallelproj projector" OFF) OPTION(DOWNLOAD_ZENODO_TEST_DATA "download zenodo data for tests" OFF) +option(DISABLE_PETSIRD "disable use of PETSIRD filetypes" OFF) option(DISABLE_UPENN "disable use of UPENN filetypes" OFF) find_package(Git QUIET) @@ -278,6 +279,41 @@ else() message(STATUS "Parallelproj projector support disabled or not available.") endif() +if(NOT DISABLE_PETSIRD) + find_package(PETSIRD CONFIG) + + if(PETSIRD_FOUND) + message(STATUS "PETSIRD support enabled (found PETSIRD).") + set(HAVE_PETSIRD ON) + + # Inspect PETSIRD transitive dependencies + if(TARGET PETSIRD::petsird) + get_target_property(_petsird_libs PETSIRD::petsird INTERFACE_LINK_LIBRARIES) + + get_target_property(_petsird_features PETSIRD::petsird INTERFACE_COMPILE_FEATURES) + get_target_property(_petsird_libs PETSIRD::petsird INTERFACE_LINK_LIBRARIES) + message(STATUS "PETSIRD::petsird dependencies: ${_petsird_libs}") + message(STATUS "PETSIRD::petsird compile features: ${_petsird_features}") + + # I have to give credit to AI for the following line. + # Probably this will not be needed, but I will keep it here for future reference. + # string(REGEX MATCH "cxx_std_([0-9]+)" _petsird_std_match "${_petsird_features}") + # if(CMAKE_MATCH_1) + # UseCXX(${CMAKE_MATCH_1}) + # endif() + + else() + message(FATAL_ERROR "PETSIRD::petsird target not found after find_package(PETSIRD)") + endif() + + else() + message(STATUS "PETSIRD not found: PETSIRD support will be disabled.") + set(HAVE_PETSIRD OFF) + endif() +else() + set(HAVE_PETSIRD OFF) +endif() + #### enable support for ctest ENABLE_TESTING() diff --git a/build_environment.yml b/build_environment.yml new file mode 100644 index 0000000000..b078dcca1a --- /dev/null +++ b/build_environment.yml @@ -0,0 +1,27 @@ +name: stir +channels: +- conda-forge +dependencies: + - cmake + - ccache + - ninja + - swig + - pre-commit>=4.0.1 + - nlohmann_json>=3.11.2 + - libboost-devel # really only libboost-headers needed, but CMake can then find a non-CONDA boost + - hdf5 # not win and not (osx and arm64) + - libitk-devel +# - cuda-version ==${{ cuda_version }} +# - cuda-driver-dev: linux +# - cuda-cudart-dev + - libparallelproj=2 + - root_base=6.34 # not win +# - zlib # win + - xtensor>=0.24.2 # PETSIRD + - xtensor-blas # PETSIRD + - howardhinnant_date>=3.0.1 # PETSIRD + - python>=3.11.3 # 3.11.3 required for PETSIRD + - numpy>=1.24.3 + - pip + - pytest + diff --git a/documentation/devel/README.md b/documentation/devel/README.md index ad8e9fb16b..3f0755ca8e 100644 --- a/documentation/devel/README.md +++ b/documentation/devel/README.md @@ -7,3 +7,4 @@ You want to contribute to STIR? Great! - Read our Git [guidelines and recommendations](git-recommendations.md) - Install git hooks (required for serious development): [instructions](git-hooks.md) - Set your editor settings appropriately: [instructions](editor-settings.md) +- To run the CI workflow locally with `act`: [instructions](local-CI.md) diff --git a/documentation/devel/local-CI.md b/documentation/devel/local-CI.md new file mode 100644 index 0000000000..7ea0b7daab --- /dev/null +++ b/documentation/devel/local-CI.md @@ -0,0 +1,37 @@ +# Running the GitHub Actions workflow locally with `act` + +The GitHub Actions CI workflow can be tested locally using [`act`](https://github.com/nektos/act). +This is useful for debugging workflow changes before pushing to GitHub. + +## Prerequisities + +Install `act` and make sure Docker or Podman is available and running. + +On Linux, the workflow can be run with an Ubuntu 24.04 + container image compatible with GitHub Actions: + +```bash +act -W .github/workflows/build-test.yml \ + -P ubuntu-24.04=ghcr.io/catthehacker/ubuntu:act-24.04 \ + --env ACT=true +``` + +* The ```-W``` option selects the workflow file to run +* The ```-P``` option maps the GitHub Actions +runner label ```ubuntu-24.04``` to a local container image. +The image in the command above is suggested online. +* The `--env ACT=true` option sets an environment variable + used by the workflow to detect that it is running under `act`. + Some GitHub steps are skipped. + +## NOTES + +* It is highly recommended to run only one job at the time. +STIR has an array of different OSes and options. +Don't try to spin them up all together in your local workstation. +* Runnning the workflow locally with ```act``` is not always identical to GitHub. +* If Docker runs out of disk space, remove old images and containers before running + +```bash +docker system prune +``` diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index aecf58f1e5..f416e29aca 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,3 +1,4 @@ +# Copyright 2025, 2026, University Medical Center Groningen # Copyright 2011-01-01 - 2011-06-30 Hammersmith Imanet Ltd # Copyright 2011-07-01 - 2012 Kris Thielemans # Copyright 2016 ETH Zurich @@ -278,6 +279,10 @@ else() target_link_libraries(stir_registries PUBLIC CUDA::cudart) endif() +if (HAVE_PETSIRD) + target_link_libraries(stir_registries PUBLIC PETSIRD::petsird) +endif() + # go and look for CMakeLists.txt files in all those directories foreach(STIR_DIR ${STIR_DIRS} ${STIR_TEST_DIRS}) ADD_SUBDIRECTORY(${STIR_DIR}) diff --git a/src/IO/CMakeLists.txt b/src/IO/CMakeLists.txt index 691258d54a..e78d97d363 100644 --- a/src/IO/CMakeLists.txt +++ b/src/IO/CMakeLists.txt @@ -80,6 +80,12 @@ if (HAVE_HDF5) ) endif() +if (HAVE_PETSIRD) +list(APPEND ${dir_LIB_SOURCES} + PETSIRDCListmodeInputFileFormat.cxx +) +endif() + endif() # MINI_STIR target_sources(${STIR_BUILDBLOCK_LIB} PRIVATE ${${dir_LIB_SOURCES}}) @@ -137,3 +143,7 @@ if (HAVE_JSON) get_target_property(TMP nlohmann_json::nlohmann_json INTERFACE_INCLUDE_DIRECTORIES) target_include_directories(${TARGET} PRIVATE "${TMP}") endif() + +if(HAVE_PETSIRD) + target_link_libraries(${TARGET} PUBLIC PETSIRD::petsird) +endif() diff --git a/src/IO/IO_registries.cxx b/src/IO/IO_registries.cxx index 242a42113c..0667f8a022 100644 --- a/src/IO/IO_registries.cxx +++ b/src/IO/IO_registries.cxx @@ -1,4 +1,5 @@ /* + Copyright 2025, 2026, University Medical Center Groningen Copyright (C) 2002-2011, Hammersmith Imanet Ltd Copyright (C) 2012, Kris Thielemans Copyright (C) 2013, Institute for Bioengineering of Catalonia @@ -66,6 +67,10 @@ # include "stir/IO/InputStreamFromROOTFileForECATPET.h" # endif +# ifdef HAVE_PETSIRD +# include "stir/IO/PETSIRDCListmodeInputFileFormat.h" +# endif + # ifdef HAVE_UPENN # include "stir/IO/PENNListmodeInputFileFormat.h" # include "stir/IO/InputStreamWithRecordsFromUPENNbin.h" @@ -156,5 +161,9 @@ static InputStreamWithRecordsFromUPENNtxt::RegisterIt dummy686062; // static RegisterInputFileFormat idummy1(2); # endif +# ifdef HAVE_PETSIRD +static RegisterInputFileFormat LMdummyPETSIRD(10); +# endif + #endif // MINI_STIR END_NAMESPACE_STIR diff --git a/src/IO/PETSIRDCListmodeInputFileFormat.cxx b/src/IO/PETSIRDCListmodeInputFileFormat.cxx new file mode 100644 index 0000000000..1374a05233 --- /dev/null +++ b/src/IO/PETSIRDCListmodeInputFileFormat.cxx @@ -0,0 +1,86 @@ +/* PETSIRDCListmodeInputFileFormat.h + + Class defining input file format for coincidence listmode data for PETSIRD. + + Copyright 2025, 2026, University Medical Center Groningen + Copyright 2025 National Physical Laboratory + + SPDX-License-Identifier: Apache-2.0 + See STIR/LICENSE.txt for details + */ +#include "stir/IO/PETSIRDCListmodeInputFileFormat.h" +#include "petsird/binary/protocols.h" +#include "petsird/hdf5/protocols.h" +#include "stir/error.h" +#include "stir/format.h" +#include +// #include "../../PETSIRD/cpp/generated/types.h" +// #include "../../PETSIRD/cpp/helpers/include/petsird_helpers.h" + +START_NAMESPACE_STIR + +/*! + + \file + \ingroup listmode + \brief Implementation of class stir::PETSIRDCListmodeInputFileFormat + + \author Nikos Efthimiou + \author Daniel Deidda + +*/ + +bool +PETSIRDCListmodeInputFileFormat::can_read(const FileSignature& signature, const std::string& filename) const +{ + + std::array hdf5_signature = { 'H', 'D', 'F', '5' }; + std::array binary_signature = { 'y', 'a', 'r', 'd' }; + + std::ifstream file(filename, std::ios::binary); + if (!file.is_open()) + { + error(format("Cannot open file: {}", filename)); + return false; + } + std::array signature_{}; + auto it = std::istreambuf_iterator(file); + auto end = std::istreambuf_iterator(); + + for (size_t i = 0; i < signature_.size() && it != end; ++i, ++it) + { + signature_[i] = *it; + } + + if (!file) + error("Stream error while reading file signature"); + + if (signature_ == hdf5_signature) + { + use_hdf5 = true; + return use_hdf5; + } + + if (signature_ == binary_signature) + { + use_hdf5 = false; + return true; + } + + // petsird::hdf5::PETSIRDReader* petsird_reader = new petsird::hdf5::PETSIRDReader(filename); + + // if (is_null_ptr(petsird_reader)) + // { + + // petsird::binary::PETSIRDReader* petsird_reader = new petsird::binary::PETSIRDReader(filename); + // if (is_null_ptr(petsird_reader)) + // { + // return false; + // } + // return true; + // } + + return false; +} + +END_NAMESPACE_STIR diff --git a/src/Testing/Temporary/CTestCostData.txt b/src/Testing/Temporary/CTestCostData.txt new file mode 100644 index 0000000000..ed97d539c0 --- /dev/null +++ b/src/Testing/Temporary/CTestCostData.txt @@ -0,0 +1 @@ +--- diff --git a/src/buildblock/CMakeLists.txt b/src/buildblock/CMakeLists.txt index 0b3bcd845f..a3106b517a 100644 --- a/src/buildblock/CMakeLists.txt +++ b/src/buildblock/CMakeLists.txt @@ -107,7 +107,12 @@ if (NOT MINI_STIR) list(APPEND ${dir_LIB_SOURCES} ProjDataGEHDF5.cxx ) -endif() + endif() + if (HAVE_PETSIRD) + list(APPEND ${dir_LIB_SOURCES} + PETSIRDInfo.cxx + ) + endif() endif() # MINI_STIR @@ -154,3 +159,7 @@ endif() if (STIR_OPENMP) target_link_libraries(${TARGET} PUBLIC ${OpenMP_EXE_LINKER_FLAGS}) endif() + +if(HAVE_PETSIRD) + target_link_libraries(${TARGET} PUBLIC PETSIRD::petsird) +endif() \ No newline at end of file diff --git a/src/buildblock/PETSIRDInfo.cxx b/src/buildblock/PETSIRDInfo.cxx new file mode 100644 index 0000000000..9ebf9e233e --- /dev/null +++ b/src/buildblock/PETSIRDInfo.cxx @@ -0,0 +1,664 @@ +/* + Copyright 2025,2026, University Medical Center Groningen + Copyright 2025, MGH / HST A. Martinos Center for Biomedical Imaging + + SPDX-License-Identifier: Apache-2.0 + See STIR/LICENSE.txt for details +*/ +/*! + +\file +\ingroup listmode +\brief implementation of class stir::CListModeDataPETSIRD + +\author Nikos Efthimiou +*/ + +#include "stir/PETSIRDInfo.h" +#include "stir/detail/PETSIRDInfo_helpers.h" + +#include "stir/warning.h" + +#include "petsird_helpers.h" +#include "petsird_helpers/create.h" +#include "petsird_helpers/geometry.h" + +#include "petsird/binary/protocols.h" +#include "petsird/hdf5/protocols.h" + +START_NAMESPACE_STIR + +/*! + \brief Infer scanner blocks and rotation axis from PETSIRD replicated modules. + \param unique_dim1_values Output set of unique translations along X. + \param unique_dim2_values Output set of unique translations along Y. + \param unique_dim3_values Output set of unique translations along Z. + \param replicated_module_list PETSIRD replicated detector modules. + \return Index of rotation axis (0=x, 1=y, 2=z) or -1 if not found. + + \details + - Extracts translation components into the provided sets. + - Uses skew-symmetric part of rotation matrices to infer axis direction. + - Emits warnings for mixed-axis rotations or inconsistencies. +*/ +int +PETSIRDInfo::figure_out_scanner_blocks_and_rotation_axis(std::set& unique_dim1_values, + std::set& unique_dim2_values, + std::set& unique_dim3_values) +{ + auto insertTranslations = [&](const petsird::RigidTransformation& trans) { + unique_dim1_values.insert(trans.matrix.at(0, 3)); + unique_dim2_values.insert(trans.matrix.at(1, 3)); + unique_dim3_values.insert(trans.matrix.at(2, 3)); + }; + + auto extractRotationMatrix = [](const petsird::RigidTransformation& trans) -> matrix::Mat3 { + matrix::Mat3 R; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + R[i][j] = trans.matrix.at(i, j); + return R; + // skew = matrix::subtract(R, matrix::transpose(R)); + // auto rot = matrix::getAxisFromSkew(skew); + }; + + std::array, 3> skew; + int detected_axis = -1; + + std::vector replicated_module_list + = petsird_scanner_info_sptr->scanner_geometry.replicated_modules; + for (const auto& module : replicated_module_list) + for (const auto& mod_trans : module.transforms) + { + insertTranslations(mod_trans); + matrix::Mat3 R = extractRotationMatrix(mod_trans); + skew = matrix::subtract(R, matrix::transpose(R)); + auto axis_vec = matrix::getAxisFromSkew(skew); + + int current_axis = -1; + for (int i = 0; i < 3; ++i) + { + if (std::abs(axis_vec[i]) > 1e-6f) + { + if (current_axis != -1) + { + warning("Rotation involves multiple axis components. Possibly non-pure rotation."); + current_axis = -2; // Sentinel for mixed axes + return -1; + } + current_axis = i; + } + } + + if (current_axis >= 0) + { + if (detected_axis == -1) + detected_axis = current_axis; + else if (detected_axis != current_axis) + warning("Inconsistent rotation axis detected between modules."); + } + } + info(format("Rotation axis of blocks inferred as axis index {}", detected_axis)); + return detected_axis; +} + +/*! + \brief Determine block element translations and radius. + \param unique_dim1_values Output set of element translations along X. + \param unique_dim2_values Output set of element translations along Y. + \param unique_dim3_values Output set of element translations along Z. + \param radius Output inferred radius (positive component orthogonal to rotation axis). + \param radius_index Output index of the axis containing the radius component. + \param rotation_axis Known rotation axis (0=x, 1=y, 2=z). + \param replicated_module_list PETSIRD replicated detector modules. + + \details + - Scans element-level transforms to infer radius and collect translations. + - Emits warnings if mixed radii are detected. +*/ +void +PETSIRDInfo::figure_out_block_element_transformations(std::set& unique_dim1_values, + std::set& unique_dim2_values, + std::set& unique_dim3_values, + float& radius, + int& radius_index, + const int rotation_axis) +{ + auto insert_translation = [&](const petsird::RigidTransformation& trans) { + unique_dim1_values.insert(trans.matrix.at(0, 3)); + unique_dim2_values.insert(trans.matrix.at(1, 3)); + unique_dim3_values.insert(trans.matrix.at(2, 3)); + }; + + auto detect_radius = [&](const petsird::RigidTransformation& trans) -> bool { + for (int i = 0; i < 3; ++i) + { + if (i == rotation_axis) + continue; + float candidate = trans.matrix.at(i, 3); + if (candidate > 0.0f) + { + radius = candidate; + radius_index = i; + return true; + } + } + return false; + }; + + for (const auto& module : petsird_scanner_info_sptr->scanner_geometry.replicated_modules) + { + for (const auto& el_trans : module.object.detecting_elements.transforms) + { + if (radius == 0.0f) + { + if (!detect_radius(el_trans)) + { + error("Unable to determine radius from translation components."); + } + } + else + { + float current = el_trans.matrix.at(radius_index, 3); + if (std::abs(current - radius) > 1e-4f) + warning("Mixed radii detected. Consider checking for misaligned modules."); + } + + insert_translation(el_trans); + } + } +} + +void +PETSIRDInfo::figure_out_block_angles(std::set& unique_angle_modules, const int rot_axis) +{ + + for (const auto& module : petsird_scanner_info_sptr->scanner_geometry.replicated_modules) + for (const auto& transform : module.transforms) + { + if (rot_axis == 0) + unique_angle_modules.insert( + std::fabs(int(1000.F * std::atan2(transform.matrix.at(1, 0), transform.matrix.at(2, 0))) / 1000.F)); + else if (rot_axis == 1) + unique_angle_modules.insert( + std::fabs(int(1000.F * std::atan2(transform.matrix.at(2, 0), transform.matrix.at(0, 0))) / 1000.F)); + else if (rot_axis == 2) + unique_angle_modules.insert( + std::fabs(int(1000.F * std::atan2(transform.matrix.at(1, 0), transform.matrix.at(0, 0))) / 1000.F)); + } +} + +PETSIRDInfo::PETSIRDInfo(const petsird::Header& header, std::string scanner_geometry) + : petsird_scanner_info_sptr(std::make_shared(header.scanner)), + petsird_header_sptr(std::make_shared(header)), + forced_geometry(scanner_geometry) +{ + + if (!petsird_scanner_info_sptr) + error("PETSIRDInfo: Null PETSIRD ScannerInformation pointer provided."); + + const auto& geom = petsird_scanner_info_sptr->scanner_geometry; + + if (geom.replicated_modules.empty()) + error("PETSIRDInfo: scanner_geometry.replicated_modules is empty."); + + if (geom.replicated_modules[0].transforms.empty()) + warning("PETSIRDInfo: replicated_modules[0].transforms is empty (rotation/angles may be unreliable)."); + + if (geom.replicated_modules[0].object.detecting_elements.transforms.empty()) + error("PETSIRDInfo: detecting_elements.transforms is empty (cannot infer element spacing/radius)."); + + //! TODO: Determine the DOI based on material + float average_doi = 0.0; + if (petsird_scanner_info_sptr->bulk_materials.size() > 0) + { + const std::string& material = petsird_scanner_info_sptr->bulk_materials[0].name; + if (material.size() > 0) + average_doi = (material == "BGO") ? 5.0f : (material == "LSO" || material == "LYSO") ? 7.0f : 0.0f; + } + + type_of_module = petsird_scanner_info_sptr->scanner_geometry.replicated_modules.size() - 1; + if (type_of_module > 0) + { + error("Multiple types of PETSIRD modules are not supported. Abort."); + } + + module_pair = petsird::TypeOfModulePair{ type_of_module, type_of_module }; + + const auto& tof_bin_edges = petsird_scanner_info_sptr->tof_bin_edges[type_of_module][type_of_module]; + info(format("Num. of TOF bins in PETSIRD {}", tof_bin_edges.NumberOfBins())); + if (tof_bin_edges.NumberOfBins() > 1) + { + info(format( + "Since the PETSIRD file has TOF information, STIR will force cylindrical geometry, as long as other things checkout.")); + forced_geometry = "cylindrical"; + } + + std::set unique_tof_values; + vector_utils::find_unique_values_2D(unique_tof_values, petsird_scanner_info_sptr->tof_resolution); + //! TODO: Supports only single type of module + numberOfModules = petsird_scanner_info_sptr->scanner_geometry.replicated_modules[0].NumberOfObjects(); + //! TODO: Supports only single type of module + numberOfElementsIndices + = petsird_scanner_info_sptr->scanner_geometry.replicated_modules[0].object.detecting_elements.NumberOfObjects(); + std::set unique_dim1_values, unique_dim2_values, unique_dim3_values; + std::set unique_tof_resolutions; + + const int rotation_axis + = figure_out_scanner_blocks_and_rotation_axis(unique_dim1_values, unique_dim2_values, unique_dim3_values); + if (rotation_axis == -1) + is_cylindrical = false; + + const std::set& main_axis = vector_utils::get_largest_vector(unique_dim1_values, unique_dim2_values, unique_dim3_values); + + int num_transaxial_blocks = numberOfModules / main_axis.size(); + info(format("I deduce that the scanner has {} transaxial number of blocks", num_transaxial_blocks)); + + float radius = 0; + int radius_indx = -1; + + std::set unique_elements_dim1_values, unique_elements_dim2_values, unique_elements_dim3_values; + figure_out_block_element_transformations( + unique_elements_dim1_values, unique_elements_dim2_values, unique_elements_dim3_values, radius, radius_indx, rotation_axis); + + std::set unique_angle_modules; + figure_out_block_angles(unique_angle_modules, rotation_axis); + + std::vector block_angular_spacing; + if (!vector_utils::get_spacing_uniform( + block_angular_spacing, unique_angle_modules, 1e-2)) /// epsilon * 10000) // relax epsilon here + { + is_cylindrical = false; + } + + std::size_t group2 = 0, group3 = 0; + { + std::vector> pet_sird_positions; + + for (uint32_t module = 0; module < 1; module++) + { + for (uint32_t elem = 0; elem < numberOfElementsIndices; elem++) + { + petsird::ExpandedDetectionBin expanded_detection_bin{ module, elem, 1 }; + auto box_shape = petsird_helpers::geometry::get_detecting_box( + *petsird_scanner_info_sptr, type_of_module, expanded_detection_bin); + CartesianCoordinate3D mean_coord; + for (auto& corner : box_shape.corners) + { // if STIR (z,y,x) -> PETSIRD (-y, -x, z) pheraps the order below needs to be changed + mean_coord.x() = +corner.c[0] / box_shape.corners.size(); + mean_coord.y() = +corner.c[1] / box_shape.corners.size(); + mean_coord.z() = +corner.c[2] / box_shape.corners.size(); + } + // mean_coord.z() += this_scanner_sptr->get_axial_crystal_spacing() / + // save mean pos into map + pet_sird_positions.push_back(mean_coord); + } + } + + inferGroupSizes_dim2_dim3(pet_sird_positions, group2, group3); + } + + std::vector element_horizontal_spacing, element_vertical_spacing; + std::set unique_elements_horizontal_values, unique_elements_vertical_values; + if (radius_indx == 0) + { + if (!vector_utils::get_spacing_uniform(element_horizontal_spacing, unique_elements_dim3_values)) + { + is_cylindrical = false; + } + if (!vector_utils::get_spacing_uniform(element_vertical_spacing, unique_elements_dim2_values)) + { + is_cylindrical = false; + } + unique_elements_horizontal_values = unique_elements_dim3_values; + unique_elements_vertical_values = unique_elements_dim2_values; + } + else + { + //! TODO: Multiple radii handling + error("Multiple radii handling not implemented yet."); + } + + { + info("Printing TOF bin edges for validation (please make sure that the STIR TOF bin edges match the PETSIRD TOF bin edges):"); + for (size_t i = 0; i < tof_bin_edges.NumberOfBins(); ++i) + { + info(format("PETSIRD TOF bin edge {}: {}", i, tof_bin_edges.edges[i]), 2); + } + } + + bool has_tile_structure = group2 > 1 && group3 > 1 && (numberOfElementsIndices % (group2 * group3) == 0); + info(format("Has tile structure: {}", has_tile_structure ? "yes" : "no")); + + if (has_tile_structure) + { + // GATE-style tiled PETSIRD + blocks_per_bucket_transaxial = group2 > 1 ? unique_elements_vertical_values.size() / group2 : group2; + info(format("blocks per bucket in transaxial direction = {}", blocks_per_bucket_transaxial)); + blocks_per_bucket_axial + = group3 > 1 ? unique_elements_horizontal_values.size() / (numberOfElementsIndices / group3) : group3; + info(format("blocks per bucket in axial direction = {}", blocks_per_bucket_axial)); + num_axial_crystals_per_block = unique_elements_horizontal_values.size() / blocks_per_bucket_axial; + num_trans_crystals_per_block = unique_elements_vertical_values.size() / blocks_per_bucket_transaxial; + } + else + { + // STIR-style flattened PETSIRD + blocks_per_bucket_transaxial = 1; + blocks_per_bucket_axial = 1; + + num_trans_crystals_per_block = unique_elements_vertical_values.size(); + num_axial_crystals_per_block = unique_elements_horizontal_values.size(); + + warning("No block structure detected: falling back to flat crystal layout."); + } + + std::vector block_axial_spacing; + vector_utils::get_spacing_uniform(block_axial_spacing, main_axis); + if (block_axial_spacing.size() < 1) + { + // std::set::iterator it = unique_elements_horizontal_values.begin(); + // std::advance(it,0); + float begin = *std::next(unique_elements_horizontal_values.begin(), 0); + // std::advance(it,unique_elements_horizontal_values.size()-1); + float end = *std::next(unique_elements_horizontal_values.begin(), unique_elements_horizontal_values.size() - 1); + block_axial_spacing.push_back(std::abs(end - begin)); + } + + info(format("I counted {} axial blocks with spacing {}", unique_dim3_values.size(), block_axial_spacing[0])); + // Check if the cyrcle area is less than 5% different from the polygon + float expected_circle_area = float(M_PI * radius * radius); + float polygon_area + = float(0.5f * unique_angle_modules.size() * radius * radius * std::sin(2.f * float(M_PI) / unique_angle_modules.size())); + info(format("Circle area: {}, Polygon area: {}, pct {}", + expected_circle_area, + polygon_area, + std::abs(expected_circle_area - polygon_area) / expected_circle_area)); + // This is the NeuroLF ratio. + forced_geometry = "BlocksOnCylindrical"; + if (std::abs(expected_circle_area - polygon_area) / expected_circle_area < 0.02f || forced_geometry == "cylindrical") + { + info(format("PETSIRDInfo: The cylindrical area {} is more than 95% matching the polygon area {}. We will presume a " + "cylindrical configuration.", + expected_circle_area, + polygon_area)); + stir_scanner_sptr.reset( + new Scanner(Scanner::User_defined_scanner, + std::string("PETSIRD_defined_scanner"), + /* num dets per ring */ + (num_transaxial_blocks * unique_elements_vertical_values.size()), + unique_dim3_values.size() * unique_elements_horizontal_values.size() /* num of rings */, + /* number of non arccor bins */ + (num_transaxial_blocks * unique_elements_vertical_values.size()) / 2, + /* number of maximum arccor bins */ + (num_transaxial_blocks * unique_elements_vertical_values.size()) / 2, + /* inner ring radius */ + radius, + /* doi */ average_doi, // average_doi, + /* ring spacing */ + element_horizontal_spacing[0], //* 10.f, + // bin_size_v + radius < 140 ? element_vertical_spacing[0] / 2 : element_vertical_spacing[0], // * 10.f,, // * 10.f, + /*intrinsic_tilt_v*/ + 0.f, + /*num_axial_blocks_per_bucket_v */ + blocks_per_bucket_axial, + /*num_transaxial_blocks_per_bucket_v*/ + blocks_per_bucket_transaxial, + /*num_axial_crystals_per_block_v*/ + num_axial_crystals_per_block, + /*num_transaxial_crystals_per_block_v*/ + num_trans_crystals_per_block, + /*num_axial_crystals_per_singles_unit_v*/ + unique_elements_horizontal_values.size() / blocks_per_bucket_axial, + /*num_transaxial_crystals_per_singles_unit_v*/ + unique_elements_vertical_values.size() / blocks_per_bucket_transaxial, + /*num_detector_layers_v*/ + 1, // num_detector_layers_v + petsird_scanner_info_sptr->energy_resolution_at_511.front(), // energy_resolution_v + 511, // reference_energy_v + tof_bin_edges.NumberOfBins(), + (tof_bin_edges.edges[1] - tof_bin_edges.edges[0]) / speed_of_light_in_mm_per_ps_div2, + *unique_tof_values.begin() * 10 // non-TOF + )); + is_cylindrical = true; + is_generic_geometry = false; + is_block_configuration = false; + } + else + { + const uint32_t forced_axial_buckets = unique_dim3_values.size(); + info("PETSIRDInfo: The cylindrical area is less than 95% matching the polygon area. We will predsume a non-cylindrical " + "configuration."); + stir_scanner_sptr.reset( + new Scanner(Scanner::User_defined_scanner, + std::string("PETSIRD_defined_scanner"), + /* num dets per ring */ + (num_transaxial_blocks * unique_elements_vertical_values.size()), + unique_dim3_values.size() * unique_elements_horizontal_values.size() /* num of rings */, + /* number of non arccor bins */ + (num_transaxial_blocks * unique_elements_vertical_values.size()) / 2, + /* number of maximum arccor bins */ + (num_transaxial_blocks * unique_elements_vertical_values.size()) / 2, + /* inner ring radius */ + radius, + /* doi */ average_doi, + /* ring spacing */ + element_horizontal_spacing[0], //* 10.f, + // bin_size_v + radius < 140 ? element_vertical_spacing[0] / 2 : element_vertical_spacing[0], // * 10.f, + /*intrinsic_tilt_v*/ + 0.f, + /*num_axial_blocks_per_bucket_v */ + forced_axial_buckets, + /*num_transaxial_blocks_per_bucket_v*/ + 1, + /*num_axial_crystals_per_block_v*/ + blocks_per_bucket_axial * num_axial_crystals_per_block, + /*num_transaxial_crystals_per_block_v*/ + blocks_per_bucket_transaxial * num_trans_crystals_per_block, + /*num_axial_crystals_per_singles_unit_v*/ + unique_elements_horizontal_values.size() / blocks_per_bucket_axial, + /*num_transaxial_crystals_per_singles_unit_v*/ + unique_elements_vertical_values.size() / blocks_per_bucket_transaxial, + /*num_detector_layers_v*/ + 1, // num_detector_layers_v + petsird_scanner_info_sptr->energy_resolution_at_511.front(), // energy_resolution_v + 511, // reference_energy_v + 1, + 0.F, + 0.F, // non-TOF + "BlocksOnCylindrical", // scanner_geometry_v + (*std::next(unique_elements_horizontal_values.begin()) + - *unique_elements_horizontal_values.begin()), // axial_crystal_spacing_v + (*std::next(unique_elements_vertical_values.begin()) + - *unique_elements_vertical_values.begin()), // transaxial_crystal_spacing_v + (*std::next(unique_elements_horizontal_values.begin()) - *unique_elements_horizontal_values.begin()) + * num_axial_crystals_per_block * blocks_per_bucket_axial, // axial_block_spacing_v + (*std::next(unique_elements_vertical_values.begin()) - *unique_elements_vertical_values.begin()) + * num_trans_crystals_per_block * blocks_per_bucket_transaxial, // transaxial_block_spacing_v + "" // crystal_map_file_name_v + )); + is_cylindrical = false; + is_generic_geometry = false; + is_block_configuration = true; + } + + /// Now let's create the PETISIRD - STIR geometry mapping + petsird_to_stir = std::make_shared(); + petsird_map_sptr = std::make_shared(); + + // enum class InnerLoopDim + //{ + // Axial, + // Tangential, + // Radial + // }; + // InnerLoopDim inner_dim = InnerLoopDim::Tangential; // determined from your groupSize analysis + + // PRECOMPUTED from previous step: + std::size_t groupSize + = blocks_per_bucket_transaxial == 1 ? 1 : num_trans_crystals_per_block; // e.g. 5, or 1 if purely monotonic + // extern InnerLoopDim inner_dim; // Axial / Tangential / Radial + + // Don't need these anymore. Keeping for future reference. + // const int num_ax = blocks_per_bucket_axial * num_axial_crystals_per_block; + // const int num_tang = blocks_per_bucket_transaxial * num_trans_crystals_per_block; + + info(format("Tile size (groupSize) = {}", groupSize)); + + // const uint32_t num_rings = static_cast(stir_scanner_sptr->get_num_rings()); + // const uint32_t num_det = static_cast(stir_scanner_sptr->get_num_detectors_per_ring()); + const uint32_t axial_blocks = static_cast(stir_scanner_sptr->get_num_axial_blocks()); + const uint32_t trans_crys = static_cast(stir_scanner_sptr->get_num_transaxial_crystals_per_block()); + const uint32_t axial_crys = static_cast(stir_scanner_sptr->get_num_axial_crystals_per_block()); + // const uint32_t layers = static_cast(stir_scanner_sptr->get_num_detector_layers()); + + for (uint32_t module = 0; module < numberOfModules; ++module) + for (uint32_t elem = 0; elem < numberOfElementsIndices; ++elem) + { + petsird::ExpandedDetectionBin bin{ module, elem, 0 }; + + int tang_pos = 0, ax_pos = 0, rad_pos = 0; + + if (!has_tile_structure) + { + // -------- STIR-origin PETSIRD (flat) -------- + const uint32_t ax_mod = module % axial_blocks; + const uint32_t tang_mod = module / axial_blocks; + + const uint32_t axial_in_block = elem % axial_crys; + const uint32_t tmp = elem / axial_crys; + const uint32_t trans_in_block = tmp % trans_crys; + const uint32_t radial = tmp / trans_crys; + + tang_pos = tang_mod * trans_crys + trans_in_block; + ax_pos = ax_mod * axial_crys + axial_in_block; + rad_pos = radial; + } + else + { + // -------- GATE-origin PETSIRD (tiled) -------- + const uint32_t groupSize = num_trans_crystals_per_block; + const uint32_t tileSize = groupSize * groupSize; + + const uint32_t tile = elem / tileSize; + const uint32_t inTile = elem % tileSize; + + const uint32_t i0 = inTile % groupSize; + const uint32_t i1 = inTile / groupSize; + + const uint32_t tang_block = tile % blocks_per_bucket_transaxial; + const uint32_t axial_block = tile / blocks_per_bucket_transaxial; + + tang_pos = tang_block * groupSize + i0 + module * (num_trans_crystals_per_block * blocks_per_bucket_transaxial); + ax_pos = axial_block * groupSize + i1; + } + + (*petsird_to_stir)[bin] = DetectionPosition<>(tang_pos, ax_pos, rad_pos); + } + + // Reverse the mapping: from STIR detpos to PETSIRD mean coord + auto map = std::make_shared(); + + for (const auto& [petsird_bin, stir_pos] : (*petsird_to_stir)) + { + auto [it, inserted] = map->emplace(stir_pos, petsird_bin); + if (!inserted) + error("Non-unique STIR DetectionPosition while building reverse map"); + } + + stir_to_petsird = map; + + if (petsird_to_stir->size() != stir_to_petsird->size()) + { + info(format("PETSIRDInfo: Map size mismatch! Forward size: {0}\n Reverse size: {1}", + petsird_to_stir->size(), + stir_to_petsird->size())); + + error("Forward and reverse maps differ in size"); + } +} + +float +PETSIRDInfo::get_detection_efficiency_for_bin(const stir::DetectionPositionPair<>& dp) const +{ + const auto& detection_bin_efficiencies = petsird_scanner_info_sptr->detection_efficiencies.detection_bin_efficiencies; + + if (!detection_bin_efficiencies) + { + return 1.f; // no efficiencies available + } + + DetectionPosition<> temp_dp1; + DetectionPosition<> temp_dp2; + + if (dp.timing_pos() < 0) + { + temp_dp1 = dp.pos2(); + temp_dp2 = dp.pos1(); + } + else + { + temp_dp1 = dp.pos1(); + temp_dp2 = dp.pos2(); + } + + auto it0 = stir_to_petsird->find(temp_dp1); + if (it0 == stir_to_petsird->end()) + { + info(format("DetectionPosition pos1(): tangential {}, axial {},radial {}", + dp.pos1().tangential_coord(), + dp.pos1().axial_coord(), + dp.pos1().radial_coord())); + error("BinNormalisationFromPETSIRD: DetectionPosition not found in STIR→PETSIRD map"); + } + + auto it1 = stir_to_petsird->find(temp_dp2); + + if (it1 == stir_to_petsird->end()) + { + info(format("DetectionPosition pos2(): tangential {}, axial {}, radial {}", + dp.pos2().tangential_coord(), + dp.pos2().axial_coord(), + dp.pos2().radial_coord())); + error("BinNormalisationFromPETSIRD: DetectionPosition not found in STIR→PETSIRD map"); + } + + const auto det0 = petsird_helpers::make_detection_bin( + *petsird_scanner_info_sptr, type_of_module, it0->second); // it0->second is ExpandedDetectionBin + + const auto det1 = petsird_helpers::make_detection_bin(*petsird_scanner_info_sptr, type_of_module, it1->second); + + return petsird_helpers::get_detection_efficiency(*petsird_scanner_info_sptr.get(), module_pair, det0, det1); +} + +float +PETSIRDInfo::get_lower_energy_threshold() const +{ + if (petsird_scanner_info_sptr->event_energy_bin_edges.size() == 0) + return 0.0f; + float min_energy = std::numeric_limits::max(); + for (const auto& bin_edges : petsird_scanner_info_sptr->event_energy_bin_edges) + { + if (bin_edges.edges.front() < min_energy) + min_energy = bin_edges.edges.front(); + } + return min_energy; +} + +float +PETSIRDInfo::get_upper_energy_threshold() const +{ + if (petsird_scanner_info_sptr->event_energy_bin_edges.size() == 0) + return 0.0f; + float max_energy = std::numeric_limits::lowest(); + for (const auto& bin_edges : petsird_scanner_info_sptr->event_energy_bin_edges) + { + if (bin_edges.edges.back() > max_energy) + max_energy = bin_edges.edges.back(); + } + return max_energy; +} + +END_NAMESPACE_STIR \ No newline at end of file diff --git a/src/buildblock/ProjDataInfo.cxx b/src/buildblock/ProjDataInfo.cxx index b72804f35b..bad5e4d14c 100644 --- a/src/buildblock/ProjDataInfo.cxx +++ b/src/buildblock/ProjDataInfo.cxx @@ -7,6 +7,7 @@ Copyright (C) 2018, University of Leeds Copyright (C) 2018, 2020-2023 University College London Copyright (C) 2016-2019, University of Hull + Copyright (C) 2025, 2026, University Medical Center Groningen This file is part of STIR. SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license @@ -169,93 +170,90 @@ ProjDataInfo::set_max_tangential_pos_num(const int max_tang_poss) max_tangential_pos_num = max_tang_poss; } -//! \todo N.E: This function is very ugly and unnessesary complicated. Could be much better. void ProjDataInfo::set_tof_mash_factor(const int new_num) { - if (scanner_ptr->is_tof_ready() && new_num > 0) + const bool tof_ready = scanner_ptr->is_tof_ready(); + + // Non-TOF mode (either scanner not TOF-ready or invalid mash factor) + if (!tof_ready || new_num <= 0) { - tof_mash_factor = new_num; - if (tof_mash_factor > scanner_ptr->get_max_num_timing_poss()) - error("ProjDataInfo::set_tof_mash_factor: TOF mashing factor (" + std::to_string(tof_mash_factor) - + +") must be smaller than or equal to the scanner's number of max timing bins (" - + std::to_string(scanner_ptr->get_max_num_timing_poss()) + ")."); - -#if 0 - // KT: code disabled as buggy but currently not needed - tof_increament_in_mm = tof_delta_time_to_mm(scanner_ptr->get_size_of_timing_pos()); - min_unmashed_tof_pos_num = - (scanner_ptr->get_max_num_timing_poss())/2; - max_unmashed_tof_pos_num = min_unmashed_tof_pos_num + (scanner_ptr->get_max_num_timing_poss()) -1; - - // Upper and lower boundaries of the timing poss; - tof_bin_unmashed_boundaries_mm.grow(min_unmashed_tof_pos_num, max_unmashed_tof_pos_num); - tof_bin_unmashed_boundaries_ps.grow(min_unmashed_tof_pos_num, max_unmashed_tof_pos_num); - - // Silently intialise the unmashed TOF bins. - for (int k = min_unmashed_tof_pos_num; k <= max_unmashed_tof_pos_num; ++k ) - { - Bin bin; - bin.timing_pos_num() = k; - // if we ever re-enable this code, there is a BUG here: - // get_k relies on num_tof_bins, so this should have been set to the unmashed value from the scanner - float cur_low = get_k(bin) - get_sampling_in_k(bin)/2.f; - float cur_high = get_k(bin) + get_sampling_in_k(bin)/2.f; - - tof_bin_unmashed_boundaries_mm[k].low_lim = cur_low; - tof_bin_unmashed_boundaries_mm[k].high_lim = cur_high; - tof_bin_unmashed_boundaries_ps[k].low_lim = static_cast(mm_to_tof_delta_time(tof_bin_unmashed_boundaries_mm[k].low_lim)); - tof_bin_unmashed_boundaries_ps[k].high_lim = static_cast(mm_to_tof_delta_time(tof_bin_unmashed_boundaries_mm[k].high_lim)); + num_tof_bins = 1; + tof_mash_factor = 0; + min_tof_pos_num = 0; + max_tof_pos_num = 0; + // we assume TOF mashing factor = 0 means non-TOF and the projector + // won't use any boundary conditions + return; + } - } -#endif - // Now, initialise the mashed TOF bins. - tof_increament_in_mm = tof_delta_time_to_mm(tof_mash_factor * scanner_ptr->get_size_of_timing_pos()); + const int max_timing_poss = scanner_ptr->get_max_num_timing_poss(); - // TODO cope with even numbers! - min_tof_pos_num = -(scanner_ptr->get_max_num_timing_poss() / tof_mash_factor) / 2; - max_tof_pos_num = min_tof_pos_num + (scanner_ptr->get_max_num_timing_poss() / tof_mash_factor) - 1; + if (new_num > max_timing_poss) + { + error("ProjDataInfo::set_tof_mash_factor: TOF mashing factor (" + std::to_string(new_num) + + ") must be smaller than or equal to the scanner's number of " + "max timing bins (" + + std::to_string(max_timing_poss) + ")."); + } + if (max_timing_poss % new_num != 0) + { + // TODO we might not need this, but are not sure. Better safe than sorry... + error("ProjDataInfo::set_tof_mash_factor: TOF mashing factor (" + std::to_string(new_num) + + ") should divide the scanner's number of max timing bins (" + std::to_string(max_timing_poss) + ")."); + } - num_tof_bins = max_tof_pos_num - min_tof_pos_num + 1; + tof_mash_factor = new_num; + // Initialise mashed TOF bins + tof_increament_in_mm = tof_delta_time_to_mm(tof_mash_factor * scanner_ptr->get_size_of_timing_pos()); - // Ensure that we have a central tof bin. - if (num_tof_bins % 2 == 0) - error("ProjDataInfo: Number of TOF bins should be an odd number. Abort."); + if (max_timing_poss % tof_mash_factor != 0) + { + error("ProjDataInfo::set_tof_mash_factor: scanner's number of timing positions (" + std::to_string(max_timing_poss) + + ") must be divisible by the TOF mashing factor (" + std::to_string(tof_mash_factor) + ")."); + } - // Upper and lower boundaries of the timing poss; - tof_bin_boundaries_mm.grow(min_tof_pos_num, max_tof_pos_num); + const int num_mashed_bins = max_timing_poss / tof_mash_factor; + num_tof_bins = num_mashed_bins; - tof_bin_boundaries_ps.grow(min_tof_pos_num, max_tof_pos_num); + // Compute min/max TOF position numbers without enforcing odd count + // + // We choose a symmetric convention: + // - If num_tof_bins is odd: bins go from -N/2 ... +N/2 + // - If num_tof_bins is even: bins go from -(N/2) ... +(N/2 - 1) + // This preserves the old behavior for odd counts and gives clean indexing for even counts. + // Note that this convention needs to match what we do in get_k(bin) + min_tof_pos_num = -num_tof_bins / 2; + max_tof_pos_num = min_tof_pos_num + num_tof_bins - 1; - for (int k = min_tof_pos_num; k <= max_tof_pos_num; ++k) - { - Bin bin; - bin.timing_pos_num() = k; - - float cur_low = get_k(bin) - get_sampling_in_k(bin) / 2.f; - float cur_high = get_k(bin) + get_sampling_in_k(bin) / 2.f; - - tof_bin_boundaries_mm[k].low_lim = cur_low; - tof_bin_boundaries_mm[k].high_lim = cur_high; - tof_bin_boundaries_ps[k].low_lim = static_cast(mm_to_tof_delta_time(tof_bin_boundaries_mm[k].low_lim)); - tof_bin_boundaries_ps[k].high_lim = static_cast(mm_to_tof_delta_time(tof_bin_boundaries_mm[k].high_lim)); - // I could imagine a better printing. - info(format("Tbin {}: {} - {} mm ({} - {} ps) = {}", - k, - tof_bin_boundaries_mm[k].low_lim, - tof_bin_boundaries_mm[k].high_lim, - tof_bin_boundaries_ps[k].low_lim, - tof_bin_boundaries_ps[k].high_lim, - get_sampling_in_k(bin))); - } - } - else if ((scanner_ptr->is_tof_ready() && new_num <= 0) - || !scanner_ptr->is_tof_ready()) // Case new_num <=, will produce non-TOF data for a TOF compatible scanner + // Upper and lower boundaries of the timing positions + tof_bin_boundaries_mm.grow(min_tof_pos_num, max_tof_pos_num); + tof_bin_boundaries_ps.grow(min_tof_pos_num, max_tof_pos_num); + + for (int k = min_tof_pos_num; k <= max_tof_pos_num; ++k) { - num_tof_bins = 1; - tof_mash_factor = 0; - min_tof_pos_num = 0; - max_tof_pos_num = 0; - // we assume TOF mashing factor = 0 means non-TOF and the projecter won't use any boundary conditions + Bin bin; + bin.timing_pos_num() = k; + + const float sampling = get_sampling_in_k(bin); + const float center = get_k(bin); + + const float cur_low = center - sampling / 2.f; + const float cur_high = center + sampling / 2.f; + + tof_bin_boundaries_mm[k].low_lim = cur_low; + tof_bin_boundaries_mm[k].high_lim = cur_high; + + tof_bin_boundaries_ps[k].low_lim = static_cast(mm_to_tof_delta_time(cur_low)); + tof_bin_boundaries_ps[k].high_lim = static_cast(mm_to_tof_delta_time(cur_high)); + + info(format("Tbin {}: {} - {} mm ({} - {} ps) = {}", + k, + tof_bin_boundaries_mm[k].low_lim, + tof_bin_boundaries_mm[k].high_lim, + tof_bin_boundaries_ps[k].low_lim, + tof_bin_boundaries_ps[k].high_lim, + sampling)); } } diff --git a/src/cmake/STIRConfig.cmake.in b/src/cmake/STIRConfig.cmake.in index cc20a4526b..38a1daa023 100644 --- a/src/cmake/STIRConfig.cmake.in +++ b/src/cmake/STIRConfig.cmake.in @@ -3,6 +3,8 @@ # Author: Kris Thielemans # Author Richard Brown # Copyright 2016, 2019, 2020, 2022, 2023 University College London +# Copyright 2025, University Medical Center Groningen +# # # Licensed under the Apache License, Version 2.0 (the "License"); # you may not use this file except in compliance with the License. @@ -158,6 +160,12 @@ if(@STIR_WITH_Parallelproj_PROJECTOR@) set(STIR_WITH_Parallelproj_PROJECTOR TRUE) endif() +if(@HAVE_PETSIRD@) + find_package(PETSIRD CONFIG REQUIRED ${STIR_FIND_TYPE}) + message(STATUS "PETSIRD support in STIR enabled.") + set(HAVE_PETSIRD TRUE) +endif() + SET(STIR_WITH_EXPERIMENTAL @STIR_ENABLE_EXPERIMENTAL@) if(STIR_WITH_EXPERIMENTAL) if(${CMAKE_VERSION} VERSION_LESS "3.12.0") diff --git a/src/cmake/STIRConfig.h.in b/src/cmake/STIRConfig.h.in index efc241dae9..5055d63a31 100644 --- a/src/cmake/STIRConfig.h.in +++ b/src/cmake/STIRConfig.h.in @@ -90,6 +90,8 @@ namespace stir { #cmakedefine HAVE_SYSTEM_GETOPT +#cmakedefine HAVE_PETSIRD + #cmakedefine STIR_DEFAULT_PROJECTOR_AS_V2 #ifndef STIR_DEFAULT_PROJECTOR_AS_V2 #define USE_PMRT diff --git a/src/include/stir/ArrayFunction.h b/src/include/stir/ArrayFunction.h index 2d7b20f172..31ec91d62e 100644 --- a/src/include/stir/ArrayFunction.h +++ b/src/include/stir/ArrayFunction.h @@ -1,6 +1,7 @@ /* Copyright (C) 2000 PARAPET partners Copyright (C) 2000- 2007, Hammersmith Imanet Ltd + Copyright (C) 2026, University Medical Center Groningen This file is part of STIR. SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license @@ -59,6 +60,7 @@ #include "stir/Array.h" #include "stir/shared_ptr.h" #include "stir/ArrayFunctionObject.h" +#include START_NAMESPACE_STIR @@ -267,6 +269,21 @@ inline void transform_array_to_periodic_indices(Array& out_array template inline void transform_array_from_periodic_indices(Array& out_array, const Array& in_array); +template +inline void find_unique_values(std::set& values, InputIt begin, InputIt end); + +// inline bool +// get_spacing_uniform(std::vector& spacing, const std::set& unsorted_block_poss, double epsilon = 1e-4) +// { +// std::vector sorted_z(unsorted_block_poss.begin(), unsorted_block_poss.end()); +// for (size_t i = 1; i < sorted_z.size(); ++i) +// { +// spacing.push_back(std::abs(sorted_z[i] - sorted_z[i - 1])); +// } + +// return std::all_of(spacing.begin(), spacing.end(), [&](float s) { return std::abs(s - spacing.front()) <= epsilon; }); +// } + END_NAMESPACE_STIR #include "stir/ArrayFunction.inl" diff --git a/src/include/stir/ArrayFunction.inl b/src/include/stir/ArrayFunction.inl index 612b7d8c03..67f01e384a 100644 --- a/src/include/stir/ArrayFunction.inl +++ b/src/include/stir/ArrayFunction.inl @@ -1,6 +1,7 @@ /* Copyright (C) 2000 PARAPET partners Copyright (C) 2000- 2007, Hammersmith Imanet Ltd + Copyright (C) 2026, University Medical Center Groningen This file is part of STIR. SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license @@ -348,4 +349,12 @@ transform_array_from_periodic_indices(Array& out_array, c } while (next(index, out_array)); } +template +inline void +find_unique_values(std::set& values, InputIt begin, InputIt end) +{ + for (auto iter = begin; iter != end; ++iter) + values.insert(*iter); +} + END_NAMESPACE_STIR diff --git a/src/include/stir/IO/PETSIRDCListmodeInputFileFormat.h b/src/include/stir/IO/PETSIRDCListmodeInputFileFormat.h new file mode 100644 index 0000000000..474831416c --- /dev/null +++ b/src/include/stir/IO/PETSIRDCListmodeInputFileFormat.h @@ -0,0 +1,62 @@ +/* + Copyright 2025, 2026 University Medical Center Groningen + Copyright 2025 National Physical Laboratory + + SPDX-License-Identifier: Apache-2.0 + See STIR/LICENSE.txt for details + */ + +/*! + + \file + \ingroup listmode + \brief Class defining input file format for coincidence listmode data for PETSIRD. + + \author Nikos Efthimiou + \author Daniel Deidda + +*/ + +#ifndef __stir_IO_PETSIRDCListmodeInputFileFormat_H__ +#define __stir_IO_PETSIRDCListmodeInputFileFormat_H__ + +// #include "boost/algorithm/string.hpp" + +#include "stir/IO/InputFileFormat.h" +#include "stir/error.h" + +#include "stir/listmode/CListModeDataPETSIRD.h" + +START_NAMESPACE_STIR + +/*! \brief Class for reading PETSIRD coincidence listmode data. + */ + +class PETSIRDCListmodeInputFileFormat : public InputFileFormat +{ +public: + const std::string get_name() const override { return "PETSIRD"; } + + //! Checks in binary data file for correct signature. + bool can_read(const FileSignature& signature, const std::string& filename) const override; + +protected: + bool actual_can_read(const FileSignature& signature, std::istream& input) const override { return false; } + + mutable bool use_hdf5 = false; + +public: + unique_ptr read_from_file(std::istream& input) const override + { + error("read_from_file for PETSIRD listmode data with istream not implemented. Sorry"); + return unique_ptr(); + } + + unique_ptr read_from_file(const std::string& filename) const override + { + info("PETSIRDCListmodeInputFileFormat: read_from_file(" + filename + ")"); + return unique_ptr(new CListModeDataPETSIRD(filename, use_hdf5)); + } +}; +END_NAMESPACE_STIR +#endif diff --git a/src/include/stir/PETSIRDInfo.h b/src/include/stir/PETSIRDInfo.h new file mode 100644 index 0000000000..4594b06d34 --- /dev/null +++ b/src/include/stir/PETSIRDInfo.h @@ -0,0 +1,171 @@ +/* + Copyright 2025, University Medical Center Groningen + + SPDX-License-Identifier: Apache-2.0 + See STIR/LICENSE.txt for detail +*/ +#ifndef __stir_listmode_PETSIRDInfo_H__ +#define __stir_listmode_PETSIRDInfo_H__ +/*! + \file PETSIRDInfo.h + \ingroup listmode + \brief Declaration of class stir::PETSIRDInfo + + \author Nikos Efthimiou +*/ + +#include "stir/DetectionPosition.h" +#include "stir/DetectionPositionPair.h" +#include "petsird/protocols.h" +#include "stir/Scanner.h" +#include "stir/DetectorCoordinateMap.h" +#include +#include + +#include "petsird_helpers.h" +#include "petsird_helpers/create.h" // for make_detection_bin +#include "petsird_helpers/geometry.h" // depending on where get_detection_efficiency lives + +START_NAMESPACE_STIR + +/*! + \brief Comparator for ordering petsird::ExpandedDetectionBin in std::map. + + \details + Orders by module_index, then element_index, then energy_index. +*/ +struct ExpandedDetectionBinLess +{ + + bool operator()(const petsird::ExpandedDetectionBin& a, const petsird::ExpandedDetectionBin& b) const + { + // Adjust field names if needed (I assume: module, element, energy_bin) + if (a.module_index < b.module_index) + return true; + if (a.module_index > b.module_index) + return false; + + if (a.element_index < b.element_index) + return true; + if (a.element_index > b.element_index) + return false; + return a.energy_index < b.energy_index; + } +}; + +/*! + \brief Mapping type from petsird::ExpandedDetectionBin to stir::DetectionPosition. +*/ +using PETSIRDToSTIRDetectorIndexMap + = std::map, ExpandedDetectionBinLess>; + +using STIRToPETSIRDDetectorIndexMap = std::map, petsird::ExpandedDetectionBin>; + +/*! + \brief Class to hold PETSIRD-related information for STIR and do any necessary conversions. +*/ +class PETSIRDInfo +{ +public: + explicit PETSIRDInfo(const petsird::Header& header, std::string scanner_geometry = ""); + + // void initialize(); + + inline std::shared_ptr get_scanner_sptr() const { return stir_scanner_sptr; } + + inline shared_ptr get_petsird_scanner_info_sptr() const { return petsird_scanner_info_sptr; } + + inline shared_ptr get_petsird_to_stir_map() const { return petsird_to_stir; } + + inline shared_ptr get_stir_to_petsird_map() const { return stir_to_petsird; } + + inline shared_ptr get_petsird_map_sptr() const { return petsird_map_sptr; } + + inline bool is_generic_geometry_used() const { return is_generic_geometry; } + + inline bool is_block_configuration_used() const { return is_block_configuration; } + + inline bool is_cylindrical_configuration_used() const { return is_cylindrical; }; + + float get_detection_efficiency_for_bin(const stir::DetectionPositionPair<>& dp) const; + + float get_lower_energy_threshold() const; + + float get_upper_energy_threshold() const; + +private: + /*! + \brief Infer scanner blocks and rotation axis from PETSIRD replicated modules. + \param unique_dim1_values Output set of unique translations along X. + \param unique_dim2_values Output set of unique translations along Y. + \param unique_dim3_values Output set of unique translations along Z. + \param replicated_module_list PETSIRD replicated detector modules. + \return Index of rotation axis (0=x, 1=y, 2=z) or -1 if not found. + */ + int figure_out_scanner_blocks_and_rotation_axis(std::set& unique_dim1_values, + std::set& unique_dim2_values, + std::set& unique_dim3_values); + /*! + \brief Compute unique module rotation angles around the given axis. + \param unique_angle_modules Output set of unique angles (radians). + \param rot_axis Rotation axis index (0=x, 1=y, 2=z). + */ + void figure_out_block_angles(std::set& unique_angle_modules, const int rot_axis); + /*! + \brief Determine block element translations and radius. + \param unique_dim1_values Output set of element translations along X. + \param unique_dim2_values Output set of element translations along Y. + \param unique_dim3_values Output set of element translations along Z. + \param radius Output inferred radius. + \param radius_index Output axis index containing the radius component. + \param rotation_axis Known rotation axis (0=x, 1=y, 2=z). + \param replicated_module_list PETSIRD replicated detector modules. +*/ + void figure_out_block_element_transformations(std::set& unique_dim1_values, + std::set& unique_dim2_values, + std::set& unique_dim3_values, + float& radius, + int& radius_index, + const int rotation_axis); + + //! Scanner information as provided by PETSIRD. + shared_ptr petsird_scanner_info_sptr; + //! Active scanner instance. + shared_ptr stir_scanner_sptr; + + shared_ptr petsird_header_sptr; + //! Number of replicated modules. + uint32_t numberOfModules; + //! Number of element indices per module. + uint32_t numberOfElementsIndices; + //! Transaxial blocks per bucket (scanner metadata). + uint32_t blocks_per_bucket_transaxial; + //! Axial blocks per bucket (scanner metadata). + uint32_t blocks_per_bucket_axial; + //! Number of axial crystals per block. + uint32_t num_axial_crystals_per_block; + //! Number of transaxial crystals per block. + uint32_t num_trans_crystals_per_block; + //! True if we should be using the cylindrical geometry + bool is_cylindrical = true; + + bool is_generic_geometry = false; + + bool is_block_configuration = false; + + petsird::TypeOfModule type_of_module; + + petsird::TypeOfModulePair module_pair; + + std::string forced_geometry = ""; + //! Mapping from PETSIRD expanded bins to STIR detection positions. + shared_ptr petsird_to_stir; + + shared_ptr stir_to_petsird; + //! Mapping from STIR detection positions to PETSIRD coordinates. + shared_ptr petsird_map_sptr; +}; + +END_NAMESPACE_STIR + +#endif \ No newline at end of file diff --git a/src/include/stir/ProjDataInfoCylindricalNoArcCorr.inl b/src/include/stir/ProjDataInfoCylindricalNoArcCorr.inl index 9f3678de3f..6627799b24 100644 --- a/src/include/stir/ProjDataInfoCylindricalNoArcCorr.inl +++ b/src/include/stir/ProjDataInfoCylindricalNoArcCorr.inl @@ -12,6 +12,7 @@ */ /* + Copyright (C) 2026, University Medical Center Groningen Copyright (C) 2000- 2005, Hammersmith Imanet Ltd This file is part of STIR. @@ -124,7 +125,7 @@ ProjDataInfoCylindricalNoArcCorr::get_bin_for_det_pair( } else { - bin.timing_pos_num() = -timing_pos_num; + bin.timing_pos_num() = -timing_pos_num - (get_num_tof_poss() % 2 == 0); return get_segment_axial_pos_num_for_ring_pair(bin.segment_num(), bin.axial_pos_num(), ring_num2, ring_num1); } } diff --git a/src/include/stir/detail/PETSIRDInfo_helpers.h b/src/include/stir/detail/PETSIRDInfo_helpers.h new file mode 100644 index 0000000000..d325cf075f --- /dev/null +++ b/src/include/stir/detail/PETSIRDInfo_helpers.h @@ -0,0 +1,280 @@ +/* + Copyright 2026, University Medical Center Groningen + Copyright 2025, MGH / HST A. Martinos Center for Biomedical Imaging + + SPDX-License-Identifier: Apache-2.0 + See STIR/LICENSE.txt for details. + + */ + +/*! + \file + \ingroup buildblock + + \brief Helper functions for PETSIRD + + \author Nikos Efthimiou +*/ + +#ifndef __stir_IO_PETSIRDInfo_helpers_H__ +#define __stir_IO_PETSIRDInfo_helpers_H__ + +#include "stir/CartesianCoordinate3D.h" +#include +#include +#include "stir/format.h" +#include "stir/info.h" +#include + +/*! + \namespace matrix + \brief Lightweight 3x3 matrix and 3D vector helpers used during PETSIRD geometry analysis. + + Provides utilities to: + - transpose a 3x3 matrix, + - subtract two 3x3 matrices, + - extract a rotation axis vector from the skew-symmetric part of a matrix. + + \details + - Mat3: std::array,3> for compact fixed-size storage. + - Vec3: std::array for simple 3D vectors. + - getAxisFromSkew(): + Given S = R - R^T (skew-symmetric part of a rotation matrix R), + returns the axis proportional to: + (S_z,y - S_y,z)/2, (S_x,z - S_z,x)/2, (S_y,x - S_x,y)/2. + For a pure rotation, S encodes the axis direction. + - These helpers assume small numerical noise; thresholds are handled by callers. + - No external dependencies; intended for quick geometric inference (e.g., rotation axis detection). + + \todo + - Might worth it to use stir::Arrays here. I started a bit of this work but lots pending. See the test_petsird_info_helpers.cxx + for some of the work done so far. + - The function are moved in ArrayFunctions +*/ + +namespace matrix +{ + +using Mat3 = std::array, 3>; +using Vec3 = std::array; + +/*! + \brief Transpose a 3x3 matrix. + \param mat Input matrix. + \return Transposed matrix. +*/ +inline Mat3 +transpose(const Mat3& mat) +{ + std::array, 3> result{}; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + result[j][i] = mat[i][j]; + return result; +} + +/*! + \brief Subtract two 3x3 matrices (A - B). + \param A Left-hand matrix. + \param B Right-hand matrix. + \return Result of A - B. +*/ +inline Mat3 +subtract(const Mat3& A, const Mat3& B) +{ + std::array, 3> result{}; + for (size_t i = 0; i < 3; ++i) + for (size_t j = 0; j < 3; ++j) + result[i][j] = A[i][j] - B[i][j]; + return result; +} + +/*! + \brief Extract rotation axis from the skew-symmetric matrix S = R - R^T. + \param S Skew-symmetric matrix. + \return Axis vector proportional to the rotation axis. +*/ +inline Vec3 +getAxisFromSkew(const Mat3& S) +{ + return { + 0.5f * (S[2][1] - S[1][2]), // x + 0.5f * (S[0][2] - S[2][0]), // y + 0.5f * (S[1][0] - S[0][1]) // z + }; +} + +} // namespace matrix + +/*! + \namespace vector_utils + \brief Helpers for spacing analysis and axis inference from coordinate sets. + + \details + - \ref vector_utils::get_spacing_uniform computes successive spacings between + sorted unique values and tests if they are uniform within a tolerance. + - \ref vector_utils::getLargestVector returns the largest of three coordinate + sets and logs which axis is inferred as “axial”. +*/ +namespace vector_utils +{ +/*! + \brief Compute spacings between sorted unique values and test uniformity. + \param spacing Output vector. Appends |x[i] - x[i-1]| for i=1..N-1, where x is the sorted version of \a unsorted_block_poss. + \param unsorted_block_poss Set of unique positions (e.g., angles or translations). + \param epsilon Tolerance for uniformity check (default 1e-4). + \return True if all spacings differ from the first spacing by <= \a epsilon, false otherwise. + + \details + - If \a spacing ends up empty (e.g., input size < 2), returns true (trivially uniform). +*/ +inline bool +get_spacing_uniform(std::vector& spacing, const std::set& unsorted_block_poss, double epsilon = 1e-4) +{ + std::vector sorted_z(unsorted_block_poss.begin(), unsorted_block_poss.end()); + for (size_t i = 1; i < sorted_z.size(); ++i) + { + spacing.push_back(std::abs(sorted_z[i] - sorted_z[i - 1])); + } + + return std::all_of(spacing.begin(), spacing.end(), [&](float s) { return std::abs(s - spacing.front()) <= epsilon; }); +} + +/*! + \brief Return the largest of three sets and report inferred axial direction. + \param x Values along X. + \param y Values along Y. + \param z Values along Z. + \return Const reference to the largest set among \a x, \a y, \a z. + + \details + Logs the index of the inferred axial direction: 0 (x), 1 (y), or 2 (z). +*/ +const std::set& +get_largest_vector(const std::set& x, const std::set& y, const std::set& z) +{ + const std::set* largest = &x; + int axis = 0; + if (y.size() > largest->size()) + { + largest = &y; + axis = 1; + } + else if (z.size() > largest->size()) + { + largest = &z; + axis = 2; + } + + // stir::info(fmt::format("I believe the axial direction is the {}.", axis)); + return *largest; +} + +/*! + \brief Collect unique values from a 1D vector. + \param values Output set for unique values. + \param input Input vector. +*/ +void +find_unique_values_1D(std::set& values, const std::vector& input) +{ + for (float val : input) + { + // std::cout << val << std::endl; + values.insert(val); + } +} + +/*! + \brief Collect unique values from a 2D vector (matrix). + \param values Output set for unique values. + \param input Input 2D vector [rows][cols]. +*/ +void +find_unique_values_2D(std::set& values, const std::vector>& input) +{ + for (size_t row = 0; row < input.size(); ++row) + for (size_t col = 0; col < input[row].size(); ++col) + values.insert(input[row][col]); +} + +} // namespace vector_utils + +bool +almostEqual(double a, double b, double tol = 1e-6) +{ + return std::fabs(a - b) <= tol; +} + +// Detect groupSize along dim2 (y) and loops along dim3 (z) +// Returns true on success, false if pattern doesn't match the assumed structure. +bool +inferGroupSizes_dim2_dim3(const std::vector>& pts, + std::size_t& groupSize_dim2, + std::size_t& groupSize_dim3, + float tol = 1e-5f) +{ + const std::size_t n = pts.size(); + groupSize_dim2 = groupSize_dim3 = 0; + + if (n == 0) + return false; + + if (n == 1) + { + groupSize_dim2 = 1; + groupSize_dim3 = 1; + return true; + } + + // STIR CartesianCoordinate3D is (z, y, x) + const float x0 = pts[0].x(); + const float z0 = pts[0].z(); + + // 1) Find how many initial points keep x and z the same + std::size_t runLen = 1; + while (runLen < n && almostEqual(pts[runLen].x(), x0, tol) && almostEqual(pts[runLen].z(), z0, tol)) + { + ++runLen; + } + + groupSize_dim2 = runLen; + + // Must tile the full array + if (groupSize_dim2 == 0 || n % groupSize_dim2 != 0) + return false; + + groupSize_dim3 = n / groupSize_dim2; + + // 2) Check each block of size groupSize_dim2 has constant x,z + for (std::size_t b = 0; b < groupSize_dim3; ++b) + { + std::size_t start = b * groupSize_dim2; + float xb = pts[start].x(); + float zb = pts[start].z(); + + for (std::size_t i = 1; i < groupSize_dim2; ++i) + { + const auto& p = pts[start + i]; + if (!almostEqual(p.x(), xb, tol) || !almostEqual(p.z(), zb, tol)) + { + return false; // pattern breaks inside a block + } + } + } + + // 3) Optionally check that z actually changes between blocks + for (std::size_t b = 1; b < groupSize_dim3; ++b) + { + float z_prev = pts[(b - 1) * groupSize_dim2].z(); + float z_curr = pts[b * groupSize_dim2].z(); + if (almostEqual(z_prev, z_curr, tol)) + { + return false; // outer loop didn't move in z + } + } + + return true; +} + +#endif \ No newline at end of file diff --git a/src/include/stir/listmode/CListModeDataBasedOnCoordinateMap.h b/src/include/stir/listmode/CListModeDataBasedOnCoordinateMap.h new file mode 100644 index 0000000000..1e9497fa42 --- /dev/null +++ b/src/include/stir/listmode/CListModeDataBasedOnCoordinateMap.h @@ -0,0 +1,56 @@ +/* + Copyright 2026, University Medical Center Groningen + Copyright 2020 Positrigo AG, Zurich + + SPDX-License-Identifier: Apache-2.0 + See STIR/LICENSE.txt for details + + */ + +/*! + + \file + \ingroup listmode + \brief Declaration of class stir::CListModeDataBasedOnCoordinateMap + + \author Nikos Efthimiou + \author Jannis Fischer +*/ + +#ifndef __stir_listmode_CListModeDataBasedOnCoordinateMap_H__ +#define __stir_listmode_CListModeDataBasedOnCoordinateMap_H__ + +#include +#include + +#include "stir/listmode/CListModeData.h" +#include "stir/DetectorCoordinateMap.h" +#include "stir/shared_ptr.h" + +// #include "stir/listmode/CListRecordSAFIR.h" + +START_NAMESPACE_STIR + +class CListModeDataBasedOnCoordinateMap : public CListModeData +{ +public: + std::string get_name() const override; + + // virtual shared_ptr> get_current_lm_file() = 0; + + SavedPosition save_get_position() override = 0; + + Succeeded set_get_position(const SavedPosition& pos) override = 0; + +protected: + std::string listmode_filename; + + mutable std::vector saved_get_positions; + virtual Succeeded open_lm_file() const = 0; + + shared_ptr map; +}; + +END_NAMESPACE_STIR + +#endif diff --git a/src/include/stir/listmode/CListModeDataPETSIRD.h b/src/include/stir/listmode/CListModeDataPETSIRD.h new file mode 100644 index 0000000000..ab0bf8df03 --- /dev/null +++ b/src/include/stir/listmode/CListModeDataPETSIRD.h @@ -0,0 +1,152 @@ +/* + + Copyright 2025, 2026, University Medical Center Groningen + Copyright 2025, MGH / HST A. Martinos Center for Biomedical Imaging + + SPDX-License-Identifier: Apache-2.0 + See STIR/LICENSE.txt for details. + + */ + +#ifndef __stir_listmode_CListModeDataPETSIRD_H__ +#define __stir_listmode_CListModeDataPETSIRD_H__ + +#include "stir/listmode/CListModeDataBasedOnCoordinateMap.h" +#include "stir/listmode/CListRecord.h" +#include "stir/listmode/CListRecordPETSIRD.h" +#include "stir/Succeeded.h" +#include "petsird_helpers.h" +#include "stir/PETSIRDInfo.h" +#include "petsird/binary/protocols.h" +#include "petsird/hdf5/protocols.h" + +START_NAMESPACE_STIR + +/*! + \class CListModeDataPETSIRD + \brief Reader for PETSIRD listmode data supporting variable geometry. + \ingroup listmode + \author Nikos Efthimiou + + \par Overview + - Supports HDF5 and binary PETSIRD formats. + - Infers scanner geometry: + - Cylindrical → creates cylindrical scanner. + - Block-based → creates block-based scanner. + - Otherwise → creates generic scanner using crystal positions. + - Builds a DetectorCoordinateMap when needed and stores to disk. + \par + + Infering the scanner geometry makes a lot of assumptions about what PET is. + In particular, it assumes: + \li A PET scanner is made of rings of detectors + \li The largest axis is the axial one. + \li So far we support only a single layer. This is partly hard-coded for simplicity. (look in the code for relevant TODOs and + comments.) + \li Some of the hardcoded assumptions are in CListRecordPETSIRD as well. + + \note Exact PETSIRD format specification is defined in the PETSIRD project documentation. + \note Initially, I wanted to: + - Is close to a cylindrical geometry ? + - then yes use a cylindrical scanner that is simpler. + - Else, is it made of blocks arranged on a cylinder. + + However, now I do the following: + - Is close to cylindrical geometry? + - yes use cylindrical scanner + - Check if blocks-on-cylinder configuration, are a good match. + - yes use blocks-on-cylinder scanner + - else use generic scanner and export the map to the disk. + + If listmode reconstruction is done, the map is regenerated on-the-fly. + +*/ +class CListModeDataPETSIRD : public CListModeDataBasedOnCoordinateMap +{ +private: + //! Snapshot of the current PETSIRD list-mode reader position. + /*! + Stores enough state to resume reading from a previously saved position + in the PETSIRD stream. This includes the current prompt/delayed stream, + the event index within the currently cached event block, the time-block + index, and optionally cached block contents. + */ + struct PetsirdCursor + { + //! Whether the cursor points to the prompt-event stream. + /*! + True for prompt events, false for delayed events. + */ + bool is_prompt = true; + //! Index of the next event within the cached event block. + std::size_t event_in_block = 0; + //! Index of the time block associated with this cursor. + std::size_t time_block_index = 0; + //! Cached PETSIRD time block at this cursor position. + petsird::TimeBlock time_block; + //! Cached PETSIRD event-time block at this cursor position. + petsird::EventTimeBlock event_block; + //! Whether time_block and event_block contain valid cached data. + bool has_cached_blocks = false; + }; + +public: + /*! + \brief Construct reader. + \param listmode_filename Path to PETSIRD listmode file. + \param use_hdf5 If true, use HDF5 reader; otherwise use binary reader. + */ + CListModeDataPETSIRD(const std::string& listmode_filename, bool use_hdf5); + + virtual shared_ptr get_empty_record_sptr() const override; + + Succeeded get_next_record(CListRecord& record_of_general_type) const override; + + SavedPosition save_get_position() override; + + Succeeded reopen_and_prime(); + + Succeeded seek_to_event_block_index(std::size_t target_event_block_index) const; + + Succeeded set_get_position(const SavedPosition& pos) override; + + virtual bool has_delayeds() const override { return m_has_delayeds; } + + Succeeded reset() override; + +protected: + virtual Succeeded open_lm_file() const override; + + shared_ptr current_lm_data_ptr; + +private: + //! Whether to use the HDF5-based PETSIRD reader. + const bool use_hdf5; + //! Index of the current event within the currently loaded event block. + mutable unsigned long int curr_event_in_event_block = 0; + //! Currently loaded PETSIRD time block. + mutable petsird::TimeBlock curr_time_block; + //! Currently loaded PETSIRD event-time block. + mutable petsird::EventTimeBlock curr_event_block; + //! Prompt/delayed classification of the current event. + //! True if the current event is a prompt event, false if it is delayed. + mutable bool curr_is_prompt = true; + //! Whether the PETSIRD data contains delayed events. + mutable bool m_has_delayeds; + //! Shared PETSIRD scanner and acquisition information. + shared_ptr petsird_info_sptr; + //! Cursor used to restore the most recently saved reader position + mutable PetsirdCursor m_saved_cursor; + //! Index of the current time block in the PETSIRD stream + mutable std::size_t m_time_block_index = 0; + //! Saved reader positions indexed by SavedPosition handles. + /*! + Each entry stores a cursor state that can later be restored, allowing + random access or rollback to previously saved positions in the PETSIRD + list-mode stream. + */ + mutable std::vector m_saved_positions; +}; + +END_NAMESPACE_STIR +#endif // CLISTMODEDATAPETSIRD_H diff --git a/src/include/stir/listmode/CListModeDataSAFIR.h b/src/include/stir/listmode/CListModeDataSAFIR.h index 66ca4344e7..be6d8e8493 100644 --- a/src/include/stir/listmode/CListModeDataSAFIR.h +++ b/src/include/stir/listmode/CListModeDataSAFIR.h @@ -1,52 +1,33 @@ -/* CListModeDataSAFIR.h +/* - Coincidence LM Data Class for SAFIR: Header File - Jannis Fischer - - Copyright 2015 ETH Zurich, Institute of Particle Physics - Copyright 2020 Positrigo AG, Zurich - - Licensed under the Apache License, Version 2.0 (the "License"); - you may not use this file except in compliance with the License. - You may obtain a copy of the License at - - http://www.apache.org/licenses/LICENSE-2.0 - - Unless required by applicable law or agreed to in writing, software - distributed under the License is distributed on an "AS IS" BASIS, - WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - See the License for the specific language governing permissions and - limitations under the License. + Copyright 2015 ETH Zurich, Institute of Particle Physics + Copyright 2020 Positrigo AG, Zurich + SPDX-License-Identifier: Apache-2.0 + See STIR/LICENSE.txt for detail. */ /*! - \file - \ingroup listmode - \brief Declaration of class stir::CListModeDataSAFIR +\file +\ingroup listmode +\brief Declaration of class stir::CListModeDataSAFIR - \author Jannis Fischer +\author Jannis Fischer */ #ifndef __stir_listmode_CListModeDataSAFIR_H__ #define __stir_listmode_CListModeDataSAFIR_H__ -#include #include -#include -#include -#include "stir/listmode/CListModeData.h" +#include "stir/listmode/CListModeDataBasedOnCoordinateMap.h" #include "stir/ProjData.h" #include "stir/ProjDataInfo.h" #include "stir/listmode/CListRecord.h" #include "stir/IO/InputStreamWithRecords.h" #include "stir/shared_ptr.h" -#include "stir/listmode/CListRecordSAFIR.h" -#include "stir/DetectorCoordinateMap.h" - START_NAMESPACE_STIR /*! @@ -57,47 +38,36 @@ START_NAMESPACE_STIR coordinates. */ template -class CListModeDataSAFIR : public CListModeData +class CListModeDataSAFIR : public CListModeDataBasedOnCoordinateMap { public: /*! Constructor - \par - Takes as arguments the filenames of the coicidence listmode file, the crystal map (text) file, and the template projection data - file - */ + \par + Takes as arguments the filenames of the coicidence listmode file, the crystal map (text) file, and the template projection data + file + */ CListModeDataSAFIR(const std::string& listmode_filename, const std::string& crystal_map_filename, const std::string& template_proj_data_filename, const double lor_randomization_sigma = 0.0); + CListModeDataSAFIR(const std::string& listmode_filename, const shared_ptr& proj_data_info_sptr); - std::string get_name() const override; shared_ptr get_empty_record_sptr() const override; Succeeded get_next_record(CListRecord& record_of_general_type) const override; + + bool has_delayeds() const override { return false; } + Succeeded reset() override; - /*! - This function should save the position in input file. This is not implemented but disabled. - Returns 0 in the moement. - \todo Maybe provide real implementation? - */ SavedPosition save_get_position() override { return static_cast(current_lm_data_ptr->save_get_position()); } - Succeeded set_get_position(const SavedPosition& pos) override { return current_lm_data_ptr->set_get_position(pos); } - /*! - Returns just false in the moment. - \todo Implement this properly to check for delayed events in LM files. - */ - bool has_delayeds() const override { return false; } + Succeeded set_get_position(const SavedPosition& pos) override { return current_lm_data_ptr->set_get_position(pos); } -private: - std::string listmode_filename; +protected: + virtual Succeeded open_lm_file() const override; mutable shared_ptr> current_lm_data_ptr; - mutable std::vector saved_get_positions; - Succeeded open_lm_file() const; - shared_ptr map; }; END_NAMESPACE_STIR - -#endif +#endif // CLISTMODEDATASAFIR_H diff --git a/src/include/stir/listmode/CListRecordPETSIRD.h b/src/include/stir/listmode/CListRecordPETSIRD.h new file mode 100644 index 0000000000..2ea0ef4bbb --- /dev/null +++ b/src/include/stir/listmode/CListRecordPETSIRD.h @@ -0,0 +1,180 @@ +/* + Copyright 2025, 2026 University Medical Center Groningen + Copyright 2025 National Physical Laboratory + + SPDX-License-Identifier: Apache-2.0 + See STIR/LICENSE.txt for details. +*/ + +/*! + + \file + \ingroup listmode + \brief Declaration of class stir::CListEventPETSIRD and stir::CListRecordPETSIRD with supporting classes + + \author Nikos Efthimiou + \author Daniel Deidda +*/ + +#ifndef __stir_listmode_CListRecordPETSIRD_H__ +#define __stir_listmode_CListRecordPETSIRD_H__ + +#include "stir/listmode/CListEventScannerWithDiscreteDetectors.h" +#include "stir/listmode/CListRecord.h" +#include "stir/DetectionPositionPair.h" +#include "stir/Succeeded.h" +#include "stir/ByteOrderDefine.h" + +#include "stir/DetectorCoordinateMap.h" +#include "stir/PETSIRDInfo.h" + +START_NAMESPACE_STIR + +template +class CListEventPETSIRD : public CListEventScannerWithDiscreteDetectors +{ +public: + inline CListEventPETSIRD(shared_ptr proj_data_info_sptr, + DetectionPositionPair<>* det_pos_pair, + bool* is_prompt) + : CListEventScannerWithDiscreteDetectors(proj_data_info_sptr), + det_pos_pair_ptr(det_pos_pair), + is_prompt_ptr(is_prompt) + {} + + // inline void get_bin(Bin& bin, const ProjDataInfo& proj_data_info) const override; + + inline bool is_prompt() const override { return *this->is_prompt_ptr; } + + bool operator==(const CListEventPETSIRD& other) const + { + if (this == &other) + return true; + + return is_prompt() == other.is_prompt() && get_detection_position() == other.get_detection_position(); + } + + inline Succeeded set_prompt(const bool prompt) override + { + *this->is_prompt_ptr = prompt; + return Succeeded::yes; + } + + virtual void get_detection_position(DetectionPositionPair<>& det_pos_pair) const override + { + det_pos_pair = *this->det_pos_pair_ptr; + } + + virtual void set_detection_position(const DetectionPositionPair<>& det_pos_pair) override + { + *this->det_pos_pair_ptr = det_pos_pair; + } + +private: + DetectionPositionPair<>* det_pos_pair_ptr = nullptr; + bool* is_prompt_ptr = nullptr; +}; + +class CListTimePETSIRD : public ListTime +{ +public: + inline unsigned long get_time_in_millisecs() const { return static_cast(time); } + inline Succeeded set_time_in_millisecs(const unsigned long time_in_millisecs) + { + time = time_in_millisecs; + return Succeeded::yes; + } + bool operator==(const CListTimePETSIRD& other) const { return time == other.time; } + inline bool is_time() const { return true; } + uint32_t time; +}; + +/*! + \brief Listmode record for PETSIRD coincidence events. + + PETSIRD stores coincidence events inside petsird::EventTimeBlock objects. + The time associated with a PETSIRD event is therefore obtained from the + enclosing EventTimeBlock, rather than from a separate time marker record. + + \par + CListModeDataPETSIRD assigns this time to the record before returning it. + Consequently, a CListRecordPETSIRD represents a coincidence event and also + carries valid timing information. + + \par + This follows the STIR convention, similar to CListRecordROOT, + that a listmode record can be both an event record and a time record. + +*/ +class CListRecordPETSIRD : public CListRecord +{ +public: + CListRecordPETSIRD(shared_ptr petsird_info_sptr, shared_ptr proj_data_info_sptr) + : event_data(make_event_data(proj_data_info_sptr, this->det_pos_pair, this->is_prompt_event)), + petsird_info_sptr(std::move(petsird_info_sptr)), + proj_data_info_sptr(std::move(proj_data_info_sptr)) + {} + + //! This record also has valid timing information. + bool is_time() const override { return true; } + //! This record represents a coincidence event. + bool is_event() const override { return true; } + + CListEvent& event() override { return *event_data; } + const CListEvent& event() const override { return *event_data; } + + CListTimePETSIRD& time() override { return time_data; } + const CListTimePETSIRD& time() const override { return time_data; } + + bool operator==(const CListRecordPETSIRD& e2) const { return event_data == e2.event_data && time_data == e2.time_data; } + + Succeeded init_from_data(petsird::CoincidenceEvent& event, const bool is_prompt = true) + { + const auto scanner_info_sptr = petsird_info_sptr->get_petsird_scanner_info_sptr(); + + const auto exp_det_0 + = petsird_helpers::expand_detection_bin(*scanner_info_sptr, + 0, // TODO type_of_module, currently we only support single module types. + event.detection_bins[0]); + const auto exp_det_1 + = petsird_helpers::expand_detection_bin(*scanner_info_sptr, + 0, // TODO type_of_module, currently we only support single module types. + event.detection_bins[1]); + auto it0 = petsird_info_sptr->get_petsird_to_stir_map()->find(exp_det_1); + auto it1 = petsird_info_sptr->get_petsird_to_stir_map()->find(exp_det_0); + if (it0 == petsird_info_sptr->get_petsird_to_stir_map()->end() || it1 == petsird_info_sptr->get_petsird_to_stir_map()->end()) + { + error("get_stir_det_pos_from_PETSIRD_id: one or both PETSIRD ids not found", + exp_det_0.module_index, + exp_det_0.element_index, + exp_det_0.energy_index); + } + + // Warning: this assumes that the PETSIRD TOF bins and the STIR ProjDataInfo + // timing positions have the same binning/mashing and number of TOF bins. + // If the STIR proj_data_info uses a different TOF mashing factor or TOF range, + // this simple offset conversion is not valid. + this->det_pos_pair = DetectionPositionPair<>( + it0->second, it1->second, static_cast(event.tof_idx) + this->proj_data_info_sptr->get_min_tof_pos_num()); + + is_prompt_event = is_prompt; + return Succeeded::yes; + } + +private: + static std::unique_ptr + make_event_data(shared_ptr proj_data_info, DetectionPositionPair<>& det_pos_pair, bool& is_prompt_event); + + std::unique_ptr event_data; + CListTimePETSIRD time_data; + + shared_ptr petsird_info_sptr; + shared_ptr proj_data_info_sptr; + + bool is_prompt_event = true; + DetectionPositionPair<> det_pos_pair; +}; + +END_NAMESPACE_STIR + +#endif diff --git a/src/include/stir/recon_buildblock/BinNormalisationFromPETSIRD.h b/src/include/stir/recon_buildblock/BinNormalisationFromPETSIRD.h new file mode 100644 index 0000000000..97273bcf71 --- /dev/null +++ b/src/include/stir/recon_buildblock/BinNormalisationFromPETSIRD.h @@ -0,0 +1,84 @@ +// +// +/*! + \file + \ingroup normalisation + + \brief Declaration of class stir::BinNormalisationFromPETSIRD + + \author Nikos Efthimiou +*/ +/* + Copyright (C) 2025, 2026, University Medical Center Groningen + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +#ifndef __stir_recon_buildblock_BinNormalisationFromPETSIRD_H__ +#define __stir_recon_buildblock_BinNormalisationFromPETSIRD_H__ + +#include "stir/recon_buildblock/BinNormalisation.h" +#include "stir/recon_buildblock/BinNormalisationWithCalibration.h" +#include "stir/RegisteredParsingObject.h" +#include "stir/shared_ptr.h" +#include "stir/PETSIRDInfo.h" + +using std::string; + +START_NAMESPACE_STIR + +class BinNormalisationFromPETSIRD + : public RegisteredParsingObject +{ +private: + using base_type = BinNormalisationWithCalibration; + +public: + //! Name which will be used when parsing a BinNormalisation object + static const char* const registered_name; + + BinNormalisationFromPETSIRD(); + + BinNormalisationFromPETSIRD(const std::string& filename); + + Succeeded set_up(const shared_ptr& exam_info_sptr, const shared_ptr&) override; + + float get_uncalibrated_bin_efficiency(const Bin& bin) const override; + + inline bool with_detector_efficiencies() const { return m_with_detector_efficiencies; } + inline bool with_dead_time() const { return m_with_dead_time; } + inline bool with_geometric_factors() const { return m_with_geometric_factors; } + +private: + bool use_hdf5; + + void set_defaults() override; + + void initialise_keymap() override; + + bool post_processing() override; + + void read_norm_data(const string& filename); + + string normalisation_filename; + + //! Flag to enable/disable detector efficiency + bool m_with_detector_efficiencies; + //! Flag to enable/disable dead time correction + bool m_with_dead_time; + //! Flag to enable/disable geometric factors + bool m_with_geometric_factors; + // shared_ptr petsird_info_sptr; + shared_ptr petsird_data_sptr; + + shared_ptr scanner_info_sptr; + + shared_ptr petsird_info_sptr; +}; + +END_NAMESPACE_STIR + +#endif diff --git a/src/listmode_buildblock/CListModeDataBasedOnCoordinateMap.cxx b/src/listmode_buildblock/CListModeDataBasedOnCoordinateMap.cxx new file mode 100644 index 0000000000..3b83ad2a05 --- /dev/null +++ b/src/listmode_buildblock/CListModeDataBasedOnCoordinateMap.cxx @@ -0,0 +1,20 @@ +/* + Copyright 2015 ETH Zurich, Institute of Particle Physics + Copyright 2020 Positrigo AG, Zurich + Copyright 2021 University College London + + SPDX-License-Identifier: Apache-2.0 + See STIR/LICENSE.txt for detail.. +*/ + +#include "stir/listmode/CListModeDataBasedOnCoordinateMap.h" + +START_NAMESPACE_STIR; + +std::string +CListModeDataBasedOnCoordinateMap::get_name() const +{ + return listmode_filename; +} + +END_NAMESPACE_STIR diff --git a/src/listmode_buildblock/CListModeDataPETSIRD.cxx b/src/listmode_buildblock/CListModeDataPETSIRD.cxx new file mode 100644 index 0000000000..a41937e171 --- /dev/null +++ b/src/listmode_buildblock/CListModeDataPETSIRD.cxx @@ -0,0 +1,318 @@ +/* + Copyright 2025, University Medical Center Groningen + Copyright 2025, MGH / HST A. Martinos Center for Biomedical Imaging + + SPDX-License-Identifier: Apache-2.0 + See STIR/LICENSE.txt for detail +*/ + +#include "stir/format.h" +#include "stir/info.h" +#include "stir/error.h" + +#include "stir/listmode/CListModeDataPETSIRD.h" + +START_NAMESPACE_STIR +/*! + \file + \ingroup listmode + \brief implementation of class stir::CListModeDataPETSIRD + + \author Daniel Deidda + \author Nikos Efthimiou +*/ + +CListModeDataPETSIRD::CListModeDataPETSIRD(const std::string& listmode_filename, bool use_hdf5) + : use_hdf5(use_hdf5) +{ + CListModeDataBasedOnCoordinateMap::listmode_filename = listmode_filename; + + petsird::Header header; + if (use_hdf5) + current_lm_data_ptr.reset(new petsird::hdf5::PETSIRDReader(listmode_filename)); + else + current_lm_data_ptr.reset(new petsird::binary::PETSIRDReader(listmode_filename)); + + m_has_delayeds = header.scanner.delayed_events_are_stored; + + current_lm_data_ptr->ReadHeader(header); + + // Get the first TimeBlock + if (!current_lm_data_ptr->ReadTimeBlocks(curr_time_block)) + error("CListModeDataPETSIRD: Could not read the first TimeBlock. Abort."); + + ++m_time_block_index; + + if (std::holds_alternative(curr_time_block)) + curr_event_block = std::get(curr_time_block); + else + error("CListModeDataPETSIRD: holds_alternative not true. Abort."); + + petsird_info_sptr = std::make_shared(header); + auto stir_scanner_sptr = petsird_info_sptr->get_scanner_sptr(); + + int tof_mash_factor = 1; + this->set_proj_data_info_sptr(std::dynamic_pointer_cast( + ProjDataInfo::construct_proj_data_info(petsird_info_sptr->get_scanner_sptr(), + 1, + petsird_info_sptr->get_scanner_sptr()->get_num_rings() - 1, + petsird_info_sptr->get_scanner_sptr()->get_num_detectors_per_ring() / 2, + petsird_info_sptr->get_scanner_sptr()->get_max_num_non_arccorrected_bins(), + /* arc_correction*/ false, + tof_mash_factor) + ->create_shared_clone())); + + shared_ptr _exam_info_sptr(new ExamInfo); + // Only PET scanners supported + _exam_info_sptr->imaging_modality = ImagingModality::PT; + _exam_info_sptr->originating_system = std::string("PETSIRD_defined_scanner"); + _exam_info_sptr->set_low_energy_thres(petsird_info_sptr->get_lower_energy_threshold()); + _exam_info_sptr->set_high_energy_thres(petsird_info_sptr->get_upper_energy_threshold()); + + this->exam_info_sptr = _exam_info_sptr; +} + +Succeeded +CListModeDataPETSIRD::open_lm_file() const +{ + // current_lm_data_ptr.reset(new petsird::hdf5::PETSIRDReader(listmode_filename)); + if (!current_lm_data_ptr->ReadTimeBlocks(curr_time_block)) + return Succeeded::no; + curr_event_block = std::get(curr_time_block); + return Succeeded::yes; +} + +shared_ptr +CListModeDataPETSIRD::get_empty_record_sptr() const +{ + shared_ptr sptr(new CListRecordPETSIRD(petsird_info_sptr, get_proj_data_info_sptr())); + return sptr; +} + +Succeeded +CListModeDataPETSIRD::get_next_record(CListRecord& record_of_general_type) const +{ + auto& record = dynamic_cast(record_of_general_type); + const auto& prompt_list = curr_event_block.prompt_events.at(0).at(0); // TODO: support multiple pairs of modules. + const auto& delayed_list = m_has_delayeds ? curr_event_block.delayed_events.at(0).at(0) : prompt_list; + + const auto& event_list = curr_is_prompt ? prompt_list : delayed_list; + + if (event_list.size() == 0) + return Succeeded::no; + + auto event = event_list.at(curr_event_in_event_block); + + if (record.init_from_data(event, curr_is_prompt) == Succeeded::no + || record_of_general_type.time().set_time_in_millisecs(curr_event_block.time_interval.start) == Succeeded::no) + { + return Succeeded::no; + } + + ++curr_event_in_event_block; + + if (curr_event_in_event_block < event_list.size()) + { + return Succeeded::yes; + } + + // - Once we hit the size of the vector + curr_event_in_event_block = 0; + + if (!m_has_delayeds || curr_is_prompt) + { + if (m_has_delayeds) + { + curr_is_prompt = false; + } + else + { + if (!current_lm_data_ptr->ReadTimeBlocks(curr_time_block)) + { + current_lm_data_ptr->Close(); + return Succeeded::no; + } + ++m_time_block_index; + curr_event_block = std::get(curr_time_block); + } + } + else + { + curr_is_prompt = true; + if (!current_lm_data_ptr->ReadTimeBlocks(curr_time_block)) + { + current_lm_data_ptr->Close(); + return Succeeded::no; + } + ++m_time_block_index; + curr_event_block = std::get(curr_time_block); + } + + return Succeeded::yes; +} + +ListModeData::SavedPosition +CListModeDataPETSIRD::save_get_position() +{ + PetsirdCursor c; + c.is_prompt = curr_is_prompt; + c.event_in_block = curr_event_in_event_block; + c.time_block_index = m_time_block_index; + + // Cache current blocks so set_get_position is instant (recommended) + c.time_block = this->curr_time_block; + c.event_block = this->curr_event_block; + c.has_cached_blocks = true; + + m_saved_positions.push_back(std::move(c)); + return static_cast(m_saved_positions.size() - 1); +} + +Succeeded +CListModeDataPETSIRD::reopen_and_prime() +{ + // ensure PETSIRD state machine is satisfied + if (current_lm_data_ptr) + { + try + { + current_lm_data_ptr->Close(); + } + catch (...) + {} + } + // current_lm_data_ptr.reset(); + if (use_hdf5) + current_lm_data_ptr.reset(new petsird::hdf5::PETSIRDReader(this->listmode_filename)); + else + current_lm_data_ptr.reset(new petsird::binary::PETSIRDReader(this->listmode_filename)); + + petsird::Header header; + current_lm_data_ptr->ReadHeader(header); + // m_eof_reached = false; + curr_event_in_event_block = 0; + curr_is_prompt = true; + m_time_block_index = 0; + // read until first EventTimeBlock + while (true) + { + if (!current_lm_data_ptr->ReadTimeBlocks(this->curr_time_block)) + { + // m_eof_reached = true; + current_lm_data_ptr->Close(); + return Succeeded::no; + } + if (std::holds_alternative(this->curr_time_block)) + { + this->curr_event_block = std::get(this->curr_time_block); + return Succeeded::yes; + } + } +} + +Succeeded +CListModeDataPETSIRD::seek_to_event_block_index(std::size_t target_event_block_index) const +{ + // assumes we are primed at event_block_index = 0 + std::size_t idx = 0; + while (idx < target_event_block_index) + { + // read next until EventTimeBlock + while (true) + { + if (!current_lm_data_ptr->ReadTimeBlocks(this->curr_time_block)) + { + // m_eof_reached = true; + current_lm_data_ptr->Close(); + return Succeeded::no; + } + if (std::holds_alternative(this->curr_time_block)) + break; + } + this->curr_event_block = std::get(this->curr_time_block); + ++idx; + } + return Succeeded::yes; +} + +Succeeded +CListModeDataPETSIRD::set_get_position(const SavedPosition& pos) +{ + if (pos >= m_saved_positions.size()) + return Succeeded::no; + const auto& c = m_saved_positions[pos]; + + // If you cached the actual blocks, you STILL must ensure the reader state + // will not be used incorrectly. Easiest: reopen+seek anyway (robust), + // then overwrite curr_* with cached data. + if (reopen_and_prime() == Succeeded::no) + return Succeeded::no; + if (seek_to_event_block_index(c.time_block_index) == Succeeded::no) + return Succeeded::no; + + // restore logical cursor + curr_is_prompt = c.is_prompt; + curr_event_in_event_block = c.event_in_block; + m_time_block_index = c.time_block_index; + if (c.has_cached_blocks) + { + this->curr_time_block = c.time_block; + this->curr_event_block = c.event_block; + } + return Succeeded::yes; +} + +Succeeded +CListModeDataPETSIRD::reset() +{ + /* \todo Not sure if this is the best way to reset the reader. + It ensures we are in a clean state, but it might be slow if the file is large and/or on a slow disk. + */ + // if (current_lm_data_ptr) + // { + // try + // { + // current_lm_data_ptr->Close(); + // } + // catch (...) + // { + // // If Close throws, treat as failure (or swallow if you must) + // return Succeeded::no; + // } + // } + + if (use_hdf5) + current_lm_data_ptr.reset(new petsird::hdf5::PETSIRDReader(this->listmode_filename)); + else + current_lm_data_ptr.reset(new petsird::binary::PETSIRDReader(this->listmode_filename)); + + curr_event_in_event_block = 0; + curr_is_prompt = true; + m_time_block_index = 0; + + try + { + while (true) + { + info(format("Reading TimeBlock index {}", m_time_block_index), 2); + if (!current_lm_data_ptr->ReadTimeBlocks(this->curr_time_block)) + return Succeeded::no; + + ++m_time_block_index; + + if (std::holds_alternative(this->curr_time_block)) + { + this->curr_event_block = std::get(this->curr_time_block); + break; + } + } + } + catch (...) + { + return Succeeded::no; + } + + return Succeeded::yes; +} + +END_NAMESPACE_STIR diff --git a/src/listmode_buildblock/CListModeDataSAFIR.cxx b/src/listmode_buildblock/CListModeDataSAFIR.cxx index 96aa9a3cb3..2256c4457c 100644 --- a/src/listmode_buildblock/CListModeDataSAFIR.cxx +++ b/src/listmode_buildblock/CListModeDataSAFIR.cxx @@ -1,22 +1,10 @@ -/* CListModeDataSAFIR.cxx +/* + Copyright 2015 ETH Zurich, Institute of Particle Physics + Copyright 2020 Positrigo AG, Zurich + Copyright 2021 University College London -Coincidence LM Data Class for SAFIR: Implementation - - Copyright 2015 ETH Zurich, Institute of Particle Physics - Copyright 2020 Positrigo AG, Zurich - Copyright 2021 University College London - - Licensed under the Apache License, Version 2.0 (the "License"); - you may not use this file except in compliance with the License. - You may obtain a copy of the License at - - http://www.apache.org/licenses/LICENSE-2.0 - - Unless required by applicable law or agreed to in writing, software - distributed under the License is distributed on an "AS IS" BASIS, - WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - See the License for the specific language governing permissions and - limitations under the License. + SPDX-License-Identifier: Apache-2.0 + See STIR/LICENSE.txt for detail */ /*! @@ -30,7 +18,6 @@ Coincidence LM Data Class for SAFIR: Implementation */ #include #include -#include #include "stir/ExamInfo.h" #include "stir/Succeeded.h" @@ -54,11 +41,11 @@ CListModeDataSAFIR::CListModeDataSAFIR(const std::string& listmode const std::string& crystal_map_filename, const std::string& template_proj_data_filename, const double lor_randomization_sigma) - : listmode_filename(listmode_filename) { + CListModeDataBasedOnCoordinateMap::listmode_filename = listmode_filename; if (!crystal_map_filename.empty()) { - map = MAKE_SHARED(crystal_map_filename, lor_randomization_sigma); + this->map = MAKE_SHARED(crystal_map_filename, lor_randomization_sigma); } else { @@ -73,7 +60,7 @@ CListModeDataSAFIR::CListModeDataSAFIR(const std::string& listmode shared_ptr template_proj_data_sptr = ProjData::read_from_file(template_proj_data_filename); this->set_proj_data_info_sptr(template_proj_data_sptr->get_proj_data_info_sptr()->create_shared_clone()); - if (open_lm_file() == Succeeded::no) + if (this->open_lm_file() == Succeeded::no) { error("CListModeDataSAFIR: Could not open listmode file " + listmode_filename + "\n"); } @@ -82,8 +69,8 @@ CListModeDataSAFIR::CListModeDataSAFIR(const std::string& listmode template CListModeDataSAFIR::CListModeDataSAFIR(const std::string& listmode_filename, const shared_ptr& proj_data_info_sptr) - : listmode_filename(listmode_filename) { + CListModeDataBasedOnCoordinateMap::listmode_filename = listmode_filename; shared_ptr _exam_info_sptr(new ExamInfo); _exam_info_sptr->imaging_modality = ImagingModality::PT; this->exam_info_sptr = _exam_info_sptr; @@ -95,13 +82,6 @@ CListModeDataSAFIR::CListModeDataSAFIR(const std::string& listmode } } -template -std::string -CListModeDataSAFIR::get_name() const -{ - return listmode_filename; -} - template shared_ptr CListModeDataSAFIR::get_empty_record_sptr() const @@ -116,10 +96,8 @@ template Succeeded CListModeDataSAFIR::get_next_record(CListRecord& record_of_general_type) const { - CListRecordT& record = static_cast(record_of_general_type); - Succeeded status = current_lm_data_ptr->get_next_record(record); - // if( status == Succeeded::yes ) record.event_SAFIR().set_map_sptr(map); - return status; + auto& record = static_cast(record_of_general_type); + return current_lm_data_ptr->get_next_record(record); } template @@ -133,14 +111,14 @@ template Succeeded CListModeDataSAFIR::open_lm_file() const { - shared_ptr stream_ptr(new fstream(listmode_filename.c_str(), ios::in | ios::binary)); + shared_ptr stream_ptr(new fstream(this->listmode_filename.c_str(), ios::in | ios::binary)); if (!(*stream_ptr)) { return Succeeded::no; } - info("CListModeDataSAFIR: opening file \"" + listmode_filename + "\"", 2); + info("CListModeDataSAFIR: opening file \"" + this->listmode_filename + "\"", 2); stream_ptr->seekg((std::streamoff)32); - current_lm_data_ptr.reset( + this->current_lm_data_ptr.reset( new InputStreamWithRecords(stream_ptr, sizeof(CListTimeDataSAFIR), sizeof(CListTimeDataSAFIR), diff --git a/src/listmode_buildblock/CListRecordPETSIRD.cxx b/src/listmode_buildblock/CListRecordPETSIRD.cxx new file mode 100644 index 0000000000..cb1c33f479 --- /dev/null +++ b/src/listmode_buildblock/CListRecordPETSIRD.cxx @@ -0,0 +1,60 @@ +/* + Copyright 2026, University Medical Center Groningen + Copyright 2026, University College London + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ +/*! + \file + \ingroup listmode + \brief Implementation of classes stir::CListRecordPETSIRD + + \author Nikos Efthimiou + \author Kris Thielemans +*/ + +#include "stir/listmode/CListRecordPETSIRD.h" +#include "stir/ProjDataInfoCylindricalNoArcCorr.h" +#include "stir/ProjDataInfoBlocksOnCylindricalNoArcCorr.h" +#include "stir/ProjDataInfoGenericNoArcCorr.h" + +START_NAMESPACE_STIR + +std::unique_ptr +CListRecordPETSIRD::make_event_data(shared_ptr proj_data_info_sptr, + DetectionPositionPair<>& det_pos_pair, + bool& is_prompt_event) +{ + // construct event of type of current ProjDataInfo + // Note: currently cumbersome due to change ProjDataInfo hierarchy. + // The following is safe... + // See https://github.com/UCL/STIR/commit/79bd05694091f7b08fb0237cb34bdbeedb256a45 + if ((proj_data_info_sptr->get_scanner_ptr()->get_scanner_geometry() == "Cylindrical") + && (dynamic_cast(proj_data_info_sptr.get()) != nullptr)) + { + return std::make_unique>( + proj_data_info_sptr, &det_pos_pair, &is_prompt_event); + } + + if ((proj_data_info_sptr->get_scanner_ptr()->get_scanner_geometry() == "BlocksOnCylindrical") + && (dynamic_cast(proj_data_info_sptr.get()) != nullptr)) + { + return std::make_unique>( + proj_data_info_sptr, &det_pos_pair, &is_prompt_event); + } + + if ((proj_data_info_sptr->get_scanner_ptr()->get_scanner_geometry() == "Generic") + && (dynamic_cast(proj_data_info_sptr.get()) != nullptr)) + { + return std::make_unique>( + proj_data_info_sptr, &det_pos_pair, &is_prompt_event); + } + + error("Unsupported ProjDataInfo type in CListRecordPETSIRD::make_event_data"); + return nullptr; +} + +END_NAMESPACE_STIR diff --git a/src/listmode_buildblock/CMakeLists.txt b/src/listmode_buildblock/CMakeLists.txt index 15badd35c4..3055c56c88 100644 --- a/src/listmode_buildblock/CMakeLists.txt +++ b/src/listmode_buildblock/CMakeLists.txt @@ -14,7 +14,8 @@ set(${dir_LIB_SOURCES} LmToProjDataWithRandomRejection.cxx CListModeDataECAT8_32bit.cxx CListRecordECAT8_32bit.cxx - CListModeDataSAFIR.cxx + CListModeDataBasedOnCoordinateMap.cxx + CListModeDataSAFIR.cxx ) if (HAVE_HDF5) @@ -24,6 +25,12 @@ list(APPEND ${dir_LIB_SOURCES} ) endif() +if (HAVE_PETSIRD) +list(APPEND ${dir_LIB_SOURCES} + CListModeDataPETSIRD.cxx +) +endif() + if (HAVE_ECAT) list(APPEND ${dir_LIB_SOURCES} CListModeDataECAT.cxx @@ -51,6 +58,10 @@ if (STIR_WITH_NiftyPET_PROJECTOR) ) endif() +if(HAVE_PETSIRD) + list(APPEND ${dir_LIB_SOURCES} + CListRecordPETSIRD.cxx) +endif() #$(dir)_REGISTRY_SOURCES:= $(dir)_registries @@ -67,3 +78,7 @@ if (HAVE_HDF5) target_include_directories(${TARGET} PRIVATE ${HDF5_INCLUDE_DIRS}) endif() endif() + +if(HAVE_PETSIRD) + target_link_libraries(${TARGET} PUBLIC PETSIRD::petsird) +endif() diff --git a/src/recon_buildblock/BinNormalisationFromPETSIRD.cxx b/src/recon_buildblock/BinNormalisationFromPETSIRD.cxx new file mode 100644 index 0000000000..61a80aacbf --- /dev/null +++ b/src/recon_buildblock/BinNormalisationFromPETSIRD.cxx @@ -0,0 +1,112 @@ +/* + Copyright (C) 2025, University Medical Center Groningen + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + See STIR/LICENSE.txt for details +*/ +/*! + \file BinNormalisationFromPETSIRD.cxx + \ingroup normalisation + + \brief Implementation for class stir::BinNormalisationFromPETSIRD + + \author Nikos Efthimiou +*/ + +#include "petsird/binary/protocols.h" +#include "petsird/hdf5/protocols.h" +#include "stir/recon_buildblock/BinNormalisationFromPETSIRD.h" +#include "stir/ProjDataInfoBlocksOnCylindricalNoArcCorr.h" +#include "stir/ProjDataInfoCylindricalNoArcCorr.h" + +START_NAMESPACE_STIR + +const char* const BinNormalisationFromPETSIRD::registered_name = "From PETSIRD"; + +void +BinNormalisationFromPETSIRD::set_defaults() +{ + base_type::set_defaults(); + normalisation_filename = ""; + use_hdf5 = false; + m_with_detector_efficiencies = true; + m_with_dead_time = true; + m_with_geometric_factors = true; +} + +void +BinNormalisationFromPETSIRD::initialise_keymap() +{ + base_type::initialise_keymap(); + parser.add_start_key("Bin Normalisation From PETSIRD"); + parser.add_key("normalisation_filename", &normalisation_filename); + parser.add_key("use hdf5", &use_hdf5); + parser.add_stop_key("End Bin Normalisation From PETSIRD"); +} + +bool +BinNormalisationFromPETSIRD::post_processing() +{ + if (base_type::post_processing()) + return true; + read_norm_data(normalisation_filename); + return false; +} + +BinNormalisationFromPETSIRD::BinNormalisationFromPETSIRD() +{ + set_defaults(); +} + +BinNormalisationFromPETSIRD::BinNormalisationFromPETSIRD(const std::string& filename) +{ + read_norm_data(filename); +} + +float +BinNormalisationFromPETSIRD::get_uncalibrated_bin_efficiency(const Bin& bin) const +{ + + DetectionPositionPair<> dp; + + if (const auto* proj_cyl = dynamic_cast(proj_data_info_sptr.get())) + { + proj_cyl->get_det_pos_pair_for_bin(dp, bin); + } + else if (const auto* proj_blk = dynamic_cast(proj_data_info_sptr.get())) + { + proj_blk->get_det_pos_pair_for_bin(dp, bin); + } + else + { + error("BinNormalisationFromPETSIRD: ProjDataInfo is neither Cylindrical nor BlocksOnCylindrical"); + } + + return petsird_info_sptr->get_detection_efficiency_for_bin(dp); +} + +Succeeded +BinNormalisationFromPETSIRD::set_up(const shared_ptr& exam_info_sptr, + const shared_ptr& proj_data_info_ptr_v) +{ + base_type::set_up(exam_info_sptr, proj_data_info_ptr_v); + + return Succeeded::yes; +} + +void +BinNormalisationFromPETSIRD::read_norm_data(const string& filename) +{ + petsird::Header header; + if (use_hdf5) + petsird_data_sptr.reset(new petsird::hdf5::PETSIRDReader(filename)); + else + petsird_data_sptr.reset(new petsird::binary::PETSIRDReader(filename)); + + petsird_data_sptr->ReadHeader(header); + + petsird_info_sptr = std::make_shared(header); +} + +END_NAMESPACE_STIR \ No newline at end of file diff --git a/src/recon_buildblock/CMakeLists.txt b/src/recon_buildblock/CMakeLists.txt index 3ddde2de2b..a3169db0cf 100644 --- a/src/recon_buildblock/CMakeLists.txt +++ b/src/recon_buildblock/CMakeLists.txt @@ -1,3 +1,4 @@ +# Copyright 2025 - University Medical Center Groningen # Copyright 2011-01-01 - 2011-06-30 Hammersmith Imanet Ltd # Copyright 2011-07-01 - 2013 Kris Thielemans @@ -105,6 +106,12 @@ list(APPEND ${dir_LIB_SOURCES} BinNormalisationFromECAT8.cxx ) +if (HAVE_PETSIRD) + list(APPEND ${dir_LIB_SOURCES} + BinNormalisationFromPETSIRD.cxx + ) +endif() + if (HAVE_HDF5) list(APPEND ${dir_LIB_SOURCES} BinNormalisationFromGEHDF5.cxx @@ -202,3 +209,7 @@ target_include_directories(${TARGET} PUBLIC #${CUVEC_INCLUDE_DIR}) if (NOT STIR_WITH_CUDA) target_compile_definitions(${TARGET} PUBLIC CUVEC_DISABLE_CUDA) endif() + +if(HAVE_PETSIRD) + target_link_libraries(${TARGET} PUBLIC PETSIRD::petsird) +endif() diff --git a/src/recon_buildblock/recon_buildblock_registries.cxx b/src/recon_buildblock/recon_buildblock_registries.cxx index c6c35a6b4d..740d0ca70e 100644 --- a/src/recon_buildblock/recon_buildblock_registries.cxx +++ b/src/recon_buildblock/recon_buildblock_registries.cxx @@ -83,6 +83,10 @@ # include "stir/recon_buildblock/CUDA/CudaRelativeDifferencePrior.h" #endif +#ifdef HAVE_PETSIRD +# include "stir/recon_buildblock/BinNormalisationFromPETSIRD.h" +#endif + #ifdef STIR_WITH_Parallelproj_PROJECTOR # include "stir/recon_buildblock/Parallelproj_projector/ForwardProjectorByBinParallelproj.h" # include "stir/recon_buildblock/Parallelproj_projector/BackProjectorByBinParallelproj.h" @@ -175,6 +179,10 @@ END_NAMESPACE_ECAT static GE::RDF_HDF5::BinNormalisationFromGEHDF5::RegisterIt dummy104; #endif +#ifdef HAVE_PETSIRD +static BinNormalisationFromPETSIRD::RegisterIt dummy105; +#endif + static FourierRebinning::RegisterIt dummyFORE; END_NAMESPACE_STIR diff --git a/src/recon_test/test_DataSymmetriesForBins_PET_CartesianGrid.cxx b/src/recon_test/test_DataSymmetriesForBins_PET_CartesianGrid.cxx index 0f68e53fa2..b31e03cb5b 100644 --- a/src/recon_test/test_DataSymmetriesForBins_PET_CartesianGrid.cxx +++ b/src/recon_test/test_DataSymmetriesForBins_PET_CartesianGrid.cxx @@ -768,7 +768,7 @@ DataSymmetriesForBins_PET_CartesianGridTests::run_tests() /*num_views=*/scanner_sptr->get_num_detectors_per_ring() / 8, /*num_tang_poss=*/64, /*arc_corrected*/ false, - /*tof_mashing*/ 116)); + /*tof_mashing*/ 117)); run_tests_for_1_projdata(proj_data_info_sptr); } @@ -783,7 +783,7 @@ DataSymmetriesForBins_PET_CartesianGridTests::run_tests() /*num_views=*/scanner_sptr->get_num_detectors_per_ring() / 8, /*num_tang_poss=*/16, /*arc_corrected*/ false, - /*tof_mashing*/ 112)); + /*tof_mashing*/ 82)); run_tests_for_1_projdata(proj_data_info_sptr); } } diff --git a/src/test/CMakeLists.txt b/src/test/CMakeLists.txt index d1d69a91e8..084f28b638 100644 --- a/src/test/CMakeLists.txt +++ b/src/test/CMakeLists.txt @@ -77,6 +77,13 @@ list (APPEND buildblock_simple_tests test_interpolate_projdata.cxx ) +if (HAVE_PETSIRD) + list (APPEND buildblock_simple_tests + test_PETSIRDInfo_helpers.cxx + ) + endif() + + endif() # MINI_STIR include(stir_test_exe_targets) diff --git a/src/test/test_ArcCorrection.cxx b/src/test/test_ArcCorrection.cxx index cb43fd3813..c4099e9406 100644 --- a/src/test/test_ArcCorrection.cxx +++ b/src/test/test_ArcCorrection.cxx @@ -176,7 +176,7 @@ ArcCorrectionTests::run_tests_tof() /*views*/ 112, /*tang_pos*/ 357, /*arc_corrected*/ false, - /*tof_mashing_factor*/ 116)); + /*tof_mashing_factor*/ 117)); cerr << "Using default range and bin-size\n"; { diff --git a/src/test/test_PETSIRDInfo_helpers.cxx b/src/test/test_PETSIRDInfo_helpers.cxx new file mode 100644 index 0000000000..4b7ec72891 --- /dev/null +++ b/src/test/test_PETSIRDInfo_helpers.cxx @@ -0,0 +1,181 @@ + +/* + Copyright (C) 2025, University Medical Center Groningen + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + See STIR/LICENSE.txt for details +*/ +/*! + + \file + \ingroup test + + \brief Test program for stir::PETSIRD hierarchy + + \author Nikos Efthimiou + + \todo The helper code is being transitioned from std::vector-based arrays + in PETSIRDInfo_helpers.h to stir::Array. This test currently checks + find_unique_values() with stir::Array inputs, while the production + PETSIRDInfo_helpers code may still use the older std::vector-based + data structures. +*/ +#include "stir/detail/PETSIRDInfo_helpers.h" +#include "stir/RunTests.h" +#include "stir/Succeeded.h" +#include "stir/Array.h" +#include "stir/make_array.h" +#include "stir/ArrayFunction.h" + +START_NAMESPACE_STIR + +class PETSIRDTests : public RunTests +{ +public: + void run_tests() override; + +private: + void test_find_unique_values_1D(); + void test_find_unique_values_2D(); + void test_get_largest_vector(); + void test_get_spacing_uniform(); + void test_get_AxisFromSkewMatrix(); + void test_infer_group_sizes_dim2_dim3(); +}; + +void +PETSIRDTests::run_tests() +{ + test_find_unique_values_1D(); + test_find_unique_values_2D(); + test_get_largest_vector(); + test_get_spacing_uniform(); + test_get_AxisFromSkewMatrix(); + test_infer_group_sizes_dim2_dim3(); +} + +void +PETSIRDTests::test_find_unique_values_1D() +{ + stir::Array<1, float> input = make_1d_array(1.0f, 2.0f, 3.0f, 2.0f, 4.0f, 1.0f, 5.0f); + stir::Array<1, float> expected = make_1d_array(1.0f, 2.0f, 3.0f, 4.0f, 5.0f); + std::set result; + + find_unique_values(result, input.begin_all_const(), input.end_all_const()); + for (const auto& val : expected) + { + this->check(result.find(val) != result.end(), format("Value {} should be in the unique set", val)); + } +} + +void +PETSIRDTests::test_find_unique_values_2D() +{ + stir::Array<2, float> input + = make_array(make_1d_array(1.0f, 2.0f, 3.0f), make_1d_array(4.0f, 2.0f, 6.0f), make_1d_array(1.0f, 8.0f, 9.0f)); + stir::Array<1, float> expected = make_1d_array(1.0f, 2.0f, 3.0f, 4.0f, 6.0f, 8.0f, 9.0f); + std::set result; + + find_unique_values(result, input.begin_all_const(), input.end_all_const()); + + for (const auto& val : expected) + { + this->check(result.find(val) != result.end(), format("Value {} should be in the unique set", val)); + } +} + +void +PETSIRDTests::test_get_largest_vector() +{ + std::set x = { 1.0f, 2.0f }; + std::set y = { 1.0f, 2.0f, 3.0f }; + std::set z = { 1.0f }; + + const std::set& largest = vector_utils::get_largest_vector(x, y, z); + this->check_if_equal(largest.size(), y.size(), "Y should be the largest vector"); +} + +void +PETSIRDTests::test_get_spacing_uniform() +{ + std::set values = { 0.0f, 2.0f, 4.0f, 6.0f, 8.0f }; + std::vector spacing; + std::set spacings; + bool is_uniform = vector_utils::get_spacing_uniform(spacing, values); + + this->check(is_uniform, "Spacing should be uniform"); + vector_utils::find_unique_values_1D(spacings, spacing); + + this->check_if_equal(static_cast(spacings.size()), 1u, "There should be one unique spacing value"); + this->check_if_equal(*spacings.begin(), 2.0f, "Spacing value should be 2.0f"); +} + +void +PETSIRDTests::test_get_AxisFromSkewMatrix() +{ + float angle_rad = static_cast(M_PI) / 6.0f; // 30 degrees + { + // Create a rotation matrix around z axis + matrix::Mat3 R = { { { 1.0f, 0.0f, 0.0f }, { 0.0f, 1.0f, 0.0f }, { 0.0f, 0.0f, 1.0f } } }; + R[0][0] = std::cos(angle_rad); + R[0][1] = -std::sin(angle_rad); + R[1][0] = std::sin(angle_rad); + R[1][1] = std::cos(angle_rad); + + std::array, 3> skew = matrix::subtract(R, matrix::transpose(R)); + auto axis = matrix::getAxisFromSkew(skew); + this->check_if_equal(axis[0], 0.0f, "X component of rotation axis should be 0"); + this->check_if_equal(axis[1], 0.0f, "Y component of rotation axis should be 0"); + this->check_if_equal(axis[2], 1.0f, "Z component of rotation axis should be 1"); + } + { + // Create a rotation matrix around x axis + matrix::Mat3 R = { { { 1.0f, 0.0f, 0.0f }, { 0.0f, 1.0f, 0.0f }, { 0.0f, 0.0f, 1.0f } } }; + R[1][1] = std::cos(angle_rad); + R[1][2] = -std::sin(angle_rad); + R[2][1] = std::sin(angle_rad); + R[2][2] = std::cos(angle_rad); + + std::array, 3> skew = matrix::subtract(R, matrix::transpose(R)); + auto axis = matrix::getAxisFromSkew(skew); + this->check_if_equal(axis[0], 1.0f, "X component of rotation axis should be 1"); + this->check_if_equal(axis[1], 0.0f, "Y component of rotation axis should be 0"); + this->check_if_equal(axis[2], 0.0f, "Z component of rotation axis should be 0"); + } +} + +void +PETSIRDTests::test_infer_group_sizes_dim2_dim3() +{ + std::vector> pts; + // Create a grid of points with groupSize_dim2 = 3 and groupSize_dim3 = 4 + for (int z = 0; z < 4; ++z) + { + for (int y = 0; y < 3; ++y) + { + // Test with float numbers and something in the x coordintate so that we don't have all zeros + pts.emplace_back(static_cast(z + 0.4), static_cast(z + 0.35), 0.1f + z / 2.f); + } + } + + std::size_t groupSize_dim2 = 0; + std::size_t groupSize_dim3 = 0; + bool success = inferGroupSizes_dim2_dim3(pts, groupSize_dim2, groupSize_dim3); + + this->check(success, "inferGroupSizes_dim2_dim3 should succeed"); + this->check_if_equal(groupSize_dim2, static_cast(3), "groupSize_dim2 should be 3"); + this->check_if_equal(groupSize_dim3, static_cast(4), "groupSize_dim3 should be 4"); +} + +END_NAMESPACE_STIR + +USING_NAMESPACE_STIR + +int +main() +{ + PETSIRDTests test; + test.run_tests(); + return test.main_return_value(); +}