diff --git a/.github/workflows/build-test.yml b/.github/workflows/build-test.yml index a39f592829..0e6d09a77e 100644 --- a/.github/workflows/build-test.yml +++ b/.github/workflows/build-test.yml @@ -38,16 +38,16 @@ jobs: strategy: matrix: include: - - os: ubuntu-latest + - os: ubuntu-24.04 compiler: gcc - compiler_version: 9 + # compiler_version: 9 cuda_version: "0" BUILD_FLAGS: "-DSTIR_OPENMP=ON" BUILD_TYPE: "Release" parallelproj: "ON" ROOT: "ON" ITK: "OFF" - - os: ubuntu-latest + - os: ubuntu-24.04 compiler: clang #compiler_version: cuda_version: "0" @@ -57,7 +57,7 @@ jobs: ROOT: "OFF" # currently using ITK 5.2 which doesn't like clang 14 ITK: "OFF" - - os: ubuntu-latest + - os: ubuntu-24.04 compiler: gcc compiler_version: 10 cuda_version: "0" @@ -66,19 +66,20 @@ jobs: parallelproj: "OFF" ROOT: "OFF" ITK: "ON" - - os: ubuntu-latest + - os: ubuntu-24.04 compiler: gcc - compiler_version: 12 + compiler_version: 14 cuda_version: "0" BUILD_FLAGS: "-DSTIR_OPENMP=ON -DCMAKE_CXX_STANDARD=20" BUILD_TYPE: "RelWithDebInfo" parallelproj: "ON" ROOT: "OFF" ITK: "ON" - - os: ubuntu-latest + - os: ubuntu-24.04 compiler: gcc - compiler_version: 12 - cuda_version: "12.1.0" + # currently CUDA doesn't support gcc 14 yet + compiler_version: 13 + cuda_version: "12.6.1" BUILD_FLAGS: "-DSTIR_OPENMP=ON -DCMAKE_CXX_STANDARD=17" BUILD_TYPE: "Release" parallelproj: "ON" @@ -86,7 +87,7 @@ jobs: ITK: "ON" - os: macOS-latest compiler: gcc - compiler_version: 11 + # compiler_version: 11 cuda_version: "0" BUILD_FLAGS: "-DSTIR_OPENMP=OFF" parallelproj: "OFF" @@ -180,7 +181,7 @@ jobs: echo CXX="$CXX" >> $GITHUB_ENV - if: matrix.cuda_version != '0' - uses: Jimver/cuda-toolkit@v0.2.11 + uses: Jimver/cuda-toolkit@v0.2.19 id: cuda-toolkit with: cuda: ${{ matrix.cuda_version }} @@ -234,13 +235,12 @@ jobs: # brew install openblas # export OPENBLAS=$(brew --prefix openblas) #python -m pip install --no-cache-dir --no-binary numpy numpy # avoid the cached .whl! - python -m pip install numpy + python -m pip install numpy pytest ;; (*) - python -m pip install numpy + python -m pip install numpy pytest ;; esac - python -m pip install pytest if test "${{matrix.parallelproj}}XX" == "ONXX"; then git clone --depth 1 --branch v1.7.3 https://github.com/gschramm/parallelproj @@ -255,13 +255,23 @@ jobs: rm -rf parallelproj fi - # Install ROOT (warning: currently only valid on Ubuntu) + # Install ROOT (warning: brittle due to OS versions etc) if test "${{matrix.ROOT}}XX" == "ONXX"; then - ROOT_file=root_v6.28.12.Linux-ubuntu20-x86_64-gcc9.4.tar.gz + case ${{matrix.os}} in + (ubuntu*) + sudo apt install libtbb-dev libvdt-dev libgif-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 + ;; + (macOS*) + ROOT_file=https://root.cern/download/root_v6.34.00.macos-15.1-arm64-clang160.tar.gz + ;; + esac wget https://root.cern/download/"$ROOT_file" tar -xzvf "$ROOT_file" rm "$ROOT_file" source root/bin/thisroot.sh + echo ROOTSYS="$ROOTSYS" >> $GITHUB_ENV fi # we'll install some dependencies with shared libraries, so need to let the OS know @@ -303,7 +313,11 @@ jobs: CMAKE_INSTALL_PREFIX=${GITHUB_WORKSPACE}/install # make available to jobs below echo CMAKE_INSTALL_PREFIX="$CMAKE_INSTALL_PREFIX" >> $GITHUB_ENV - EXTRA_BUILD_FLAGS="-DBUILD_SWIG_PYTHON=ON -DPYTHON_EXECUTABLE=`which python`" + 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" + 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} -DDISABLE_STIR_LOCAL=OFF -DSTIR_LOCAL=${GITHUB_WORKSPACE}/examples/C++/using_STIR_LOCAL" @@ -321,6 +335,13 @@ jobs: ;; esac + # Enable tmate debugging of manually-triggered workflows if the input option was provided + - name: Setup tmate session if triggered + #if: ${{ failure() }} + uses: mxschmitt/action-tmate@v3 + timeout-minutes: 30 + if: ${{ github.event_name == 'workflow_dispatch' && inputs.debug_enabled == 'true' }} + - name: build shell: bash env: @@ -330,12 +351,6 @@ jobs: source ${GITHUB_WORKSPACE}/my-env/bin/activate cmake --build . -j 2 --config ${BUILD_TYPE}} --target install - # Enable tmate debugging of manually-triggered workflows if the input option was provided - - name: Setup tmate session if triggered - uses: mxschmitt/action-tmate@v3 - timeout-minutes: 15 - if: ${{ github.event_name == 'workflow_dispatch' && inputs.debug_enabled == 'true' }} - - name: ctest shell: bash env: diff --git a/.mailmap b/.mailmap index e3444addf2..9f3aaa1105 100644 --- a/.mailmap +++ b/.mailmap @@ -62,6 +62,7 @@ Georg Schramm <40211162+gschramm@users.noreply.githu Georg Schramm Georg Schramm Nicole Jurjew +Nicole Jurjew <86268422+NicoleJurjew@users.noreply.github.com> Tahereh Niknejad Tahereh Niknejad Sam D Porter <92305641+samdporter@users.noreply.github.com> @@ -70,3 +71,4 @@ Matthew Strugari Matthew Strugari <56315593+mastergari@users.noreply.github.com> Imraj Singh Imraj Singh <72553490+Imraj-Singh@users.noreply.github.com> +Dimitra Kyriakopoulou diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index d77c551dff..508b1dfee3 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -1,6 +1,6 @@ repos: -- repo: https://github.com/doublify/pre-commit-clang-format - rev: '62302476' - hooks: - - id: clang-format - files: \.(c|cc|cxx|cpp|cu|h|hpp|hxx|inl|txx)$ +- repo: https://github.com/pre-commit/mirrors-clang-format + rev: v14.0.6 + hooks: + - id: clang-format + files: \.(c|cc|cxx|cpp|cu|h|hpp|hxx|inl|txx)$ diff --git a/documentation/credits.htm b/documentation/credits.htm index c409ad0ee6..527908494c 100644 --- a/documentation/credits.htm +++ b/documentation/credits.htm @@ -142,6 +142,7 @@

contributions since the PARAPET EU project (April 2000)

Emond, Elise Falcón, Carlés Twyman, Robert + Singh, Imraj Jurjew, Nicole Bertolli, Ottavia Tsai, Yu-jung @@ -244,7 +245,7 @@

Web site

-Last modified: Wed Feb 07 2024 +Last modified: Tue Jul 23 2024 diff --git a/documentation/devel/editor-settings.md b/documentation/devel/editor-settings.md index 999f00a84d..5ebda5e4ce 100644 --- a/documentation/devel/editor-settings.md +++ b/documentation/devel/editor-settings.md @@ -1,5 +1,9 @@ # Developer documentation: editor settings White-spaces and indentation with multiple developers are a pain. Please adhere to -our white-space policy, which we try to enforce via [clang-format](https://clang.llvm.org/docs/ClangFormat.html). -Check that site for integration with your editor/IDE. +our white-space policy, which we enforce via [pre-commit](https://pre-commit.com/), including +running of [clang-format](https://clang.llvm.org/docs/ClangFormat.html). +See [git-hooks.md](git-hooks.md) for more information. + +Developer experience will be best if you set your editor/IDE to use the same settings, +including versions of the formatters (`clang-format` in particular). diff --git a/documentation/devel/git-hooks.md b/documentation/devel/git-hooks.md index 6bb563c64a..dd76c7a206 100644 --- a/documentation/devel/git-hooks.md +++ b/documentation/devel/git-hooks.md @@ -3,11 +3,27 @@ We use [pre-commit](https://pre-commit.com) to check coding conventions, including white-space formatting rules. -Unfortunately, exact formatting used depends on the `clang-format` version. Therefore, -you do have to install the same version as what is currently used in our GitHub Action. +Unfortunately, exact formatting used depends on the `clang-format` (and others?) version. Therefore, +you might need to install the same version as what is currently used in our +[.pre-commit-config.yaml](../../.pre-commit-config.yaml) for your editor to pick up +the same version. ## Installation of software +### TL;DR + +```sh +pip/conda install pre-commit +git add blabla +pre-commit run +# optionally check files +# "add" changes made by the hooks +git add blabla +git commit +``` + +### More detail + We highly recommend to use `conda`: ```sh cd /whereever/STIR @@ -48,4 +64,4 @@ or one-off with git commit --no-verify -You will then need to run `pre-commit run --all-files` when done before we will merge your pull request. \ No newline at end of file +You will then need to run `pre-commit run --all-files` when done before we will merge your pull request. diff --git a/documentation/release_6.2.htm b/documentation/release_6.2.htm index 10c7f56ccf..fa8387c9f9 100644 --- a/documentation/release_6.2.htm +++ b/documentation/release_6.2.htm @@ -27,7 +27,8 @@

Overall summary

This release contains mainly code written by Nicole Jurjew (UCL) (SSRB for TOF), - Imraj Singh (UCL) (CUDA version of the Relative Difference Prior) and Kris Thielemans (UCL). + Imraj Singh (UCL) (CUDA version of the Relative Difference Prior), + Markus Jehl (Positrigo) (fixes for blocks on cylindrical) and Kris Thielemans (UCL).

Patch release info

diff --git a/documentation/release_6.3.htm b/documentation/release_6.3.htm new file mode 100644 index 0000000000..dc6cbad024 --- /dev/null +++ b/documentation/release_6.3.htm @@ -0,0 +1,123 @@ + + + + Summary of changes in STIR release 6.3 + + + +

Summary of changes in STIR release 6.3

+ + +

Overall summary

+ + +

Patch release info

+ + +

Summary for end users (also to be read by developers)

+ + +

New functionality

+
    +
  • + The analytic Spline Reconstruction Technique (SRT) algorithm has been added in 2 different versions: for PET + (inverting the 2D Radon transform) and for SPECT (inverting the 2D attenuated Radon transform). The latter + allows quantitatively correct analytic reconstruction of SPECT data (after scatter correction). + + The reference for the implemented algorithms is:
    + Fokas, A. S., A. Iserles, and V. Marinakis. + Reconstruction algorithm for single photon emission computed tomography and its numerical implementation. + Journal of the Royal Society Interface* 3.6 (2006): 45-54. +
    + The reference for the implementation is:
    + Dimitra Kyriakopoulou, Analytical and Numerical Aspects of Tomography, UCL PhD thesis, (not yet publically accessible) +
    + PR #1420 +
    • +
  • + ScatterSimulation can now downsample the scanner transaxially (crystals per ring) for BlocksOnCylindrical, + scanners, which speeds up ScatterEstimation considerably. By default, downsampling the detectors per reading + is disabled for backwards compatibility.
    + PR #1291 +
    • +
  • + Data from GE Discovery MI systems in RDF9 should now be readable. TOF information on these scanners has also been added. + However, projection data is currently still always returned as non-TOF (but list-mode data is read as TOF).
    + PR #1503 +
    • +
  • + stir_timings has now an extra option to parse a par-file for a projector-pair. +
    • +
  • + Added the ability to set a forward projector for mask projection in the ScatterEstimation class.
    + PR #1530 +
    • +
  • + Duration in sinogram interfile/exam_info obtained from lm_to_projdata has the correct value if we unlist all the events. This is not true for ROOT files
    + PR #1519 +
    • +
+ + +

Changed functionality

+
    +
  • + Default ECAT scanner configurations updated to use a negative intrinsic tilt. +
    • +
+ +

Bug fixes

+
    +
  • + Fixed a bug in the scatter estimation code (introduced in release 5.1.0) if input data is 3D and "cylindrical" (there was no bug for "blocksoncylindrical" data). + The scatter estimation runs on data constructed via SSRB. However, the attenuation correction factors were incorrectly obtained with adding of oblique segments (as opposed to averaging). + This resulted in intermediate images that had the wrong attenuation correction which were approximately num_segments times larger. This was compensated by the tail-fitting, but resulted in unexpected scale factors (scale factors were around 1/num_segments times what was expected). + This means that if you used the "min/max scale factor" feature in the scatter estimate, you will have to adjust your threshold values. Expected scatter tail-fitting scale factors should now be restored to ~1-1.5 (depending on the amount of multiple and out-of-FOV scatter). + See Issue #1532 for more detail. Fixed by using averaging functionality of SSRB instead of adding segments for attenuation correction factors. + PR #1531 +
    • +
+ +

Build system

+
    +
  • Enable more diagnostics in CMake when finding CERN's ROOT (we used to silence them)
    + PR #1552 +
    • +
+ +

Known problems

+

See our issue tracker.

+ + +

What is new for developers (aside from what should be obvious from the above):

+ +

New functionality

+ +
  • + ProjDataInMemory read_from_file method now returns a ProjDataInMemory object. +
    • + +

    Changed functionality

    + + +

    Bug fixes

    +
      +
    • Fixed minor incompatibility with gcc-14 and clang-18 buy adding an extra include file
      + PR #1552 +
      • +
    + +

    Other code changes

    + + +

    Test changes

    + + +

    C++ tests

    + + +

    recon_test_pack

    + + + + diff --git a/examples/python/projdata_visualisation/ProjDataVisualisation.py b/examples/python/projdata_visualisation/ProjDataVisualisation.py index 983d50e62e..5d66524a8f 100644 --- a/examples/python/projdata_visualisation/ProjDataVisualisation.py +++ b/examples/python/projdata_visualisation/ProjDataVisualisation.py @@ -182,7 +182,6 @@ def update_display_image(self): f"TOF bin: {self.UI_groupbox_projdata_dimensions.value(ProjDataDims.TIMING_POS)}") ax.yaxis.set_label_text("Views/projection angle") ax.xaxis.set_label_text("Tangential positions") - ax.xaxis.set_label_text("TOF bins") elif self.viewgram_radio_button.isChecked(): image = self.get_viewgram_numpy_array() ax.title.set_text( @@ -191,7 +190,6 @@ def update_display_image(self): f"TOF bin: {self.UI_groupbox_projdata_dimensions.value(ProjDataDims.TIMING_POS)}") ax.yaxis.set_label_text("Axial positions") ax.xaxis.set_label_text("Tangential positions") - ax.xaxis.set_label_text("TOF bins") else: msg = f"Error: No radio button is checked... How did you get here?\n" raise Exception(msg) diff --git a/examples/samples/SRT2D.par b/examples/samples/SRT2D.par new file mode 100644 index 0000000000..4224bd4c01 --- /dev/null +++ b/examples/samples/SRT2D.par @@ -0,0 +1,18 @@ +SRT2DParameters := +input file := sino_hof_no_noise.hs + +xy output image size (in pixels) := -1 +z output image size (in pixels) := -1 +zoom := 1 + +;post-filter type := Wiener +;Wiener Filter Parameters := +;End Wiener Filter Parameters := + +;post-filter type := Gamma +;Gamma Filter Parameters := +;End Gamma Filter Parameters := + +output filename prefix := SRT2D + +END := diff --git a/examples/samples/SRT2DSPECT.par b/examples/samples/SRT2DSPECT.par new file mode 100644 index 0000000000..798038cbbe --- /dev/null +++ b/examples/samples/SRT2DSPECT.par @@ -0,0 +1,20 @@ +SRT2DSPECTParameters := + +;SRT2DSPECT receives as input the forward projection of the attenuation image (attensino.hs) and the measured emission sinogram (sino.hs). +input file := sino.hs +attenuation projection filename := attensino.hs + +xy output image size (in pixels) := -1 + +;post-filter type := Wiener +;Wiener Filter Parameters := +;End Wiener Filter Parameters := + +;post-filter type := Gamma +;Gamma Filter Parameters := +;End Gamma Filter Parameters := + +output filename prefix := srt_recon + + +end:= diff --git a/recon_test_pack/SPECT_test_Interfile_header.hs b/recon_test_pack/SPECT_test_Interfile_header.hs new file mode 100644 index 0000000000..4864542529 --- /dev/null +++ b/recon_test_pack/SPECT_test_Interfile_header.hs @@ -0,0 +1,62 @@ +!INTERFILE := +; This is a sample minimal header for SPECT tomographic data +; The format is as per the 3.3 Interfile standard (aside from time frame info) + +!imaging modality := nucmed + +; name of file with binary data +name of data file := SPECT_test_Interfile_header.s + +!version of keys := 3.3 +!GENERAL DATA := +!GENERAL IMAGE DATA := +!type of data := Tomographic + +; optional keywords specifying patient position +; patient rotation := prone +; patient orientation := feet_in + +imagedata byte order := LITTLEENDIAN + +number of energy windows:=1 +energy window lower level[1]:=120 +energy window upper level[1]:=160 + +!SPECT STUDY (General) := +; specify how the data are stored on disk +; here given as "single-precision float" (you could have "unsigned integer" data instead) +!number format := float +!number of bytes per pixel := 4 +!number of projections := 120 +; total rotation (or coverage) angle (in degrees) +!extent of rotation := 360 +process status := acquired +!SPECT STUDY (acquired data):= +; rotation info (e.g. clock-wise or counter-clock wise) +!direction of rotation := CW +start angle := 180 + +; Orbit definition +orbit := Circular +; radius in mm +Radius := 166.5 +; or +; orbit := Non-circular +; give a list of "radii", one for every position +; Radii := {150, 151, 153, ....} + +; pixel sizes in the acquired data, first in "transverse" direction, then in "axial" direction +!matrix size [1] := 111 +!scaling factor (mm/pixel) [1] := 3 +!matrix size [2] := 47 +!scaling factor (mm/pixel) [2] := 3.27 + +; optional keywords specifying frame duration etc +; These are not according to the Interfile 3.3 specification +; Currently only useful in STIR for dynamic applications +; (but a "time frame" is considered to be all projections acquired at the same time) +;number of time frames := 1 +;image duration (sec)[1] := 0 +;image relative start time (sec)[1] := 0 + +!END OF INTERFILE := diff --git a/recon_test_pack/SPECT_test_Interfile_header.s b/recon_test_pack/SPECT_test_Interfile_header.s new file mode 100644 index 0000000000..8b13789179 --- /dev/null +++ b/recon_test_pack/SPECT_test_Interfile_header.s @@ -0,0 +1 @@ + diff --git a/recon_test_pack/SRT2DSPECT_test_sim.par b/recon_test_pack/SRT2DSPECT_test_sim.par new file mode 100644 index 0000000000..58e35c6a5c --- /dev/null +++ b/recon_test_pack/SRT2DSPECT_test_sim.par @@ -0,0 +1,20 @@ +SRT2DSPECTParameters := +; test file for SRT2DSPECT +input file := my_sino${suffix}.hs +attenuation projection filename := my_attenuation_sino${suffix}.hs + +xy output image size (in pixels) := 91 +zoom := 0.5 + +;post-filter type := Wiener +;Wiener Filter Parameters := +;End Wiener Filter Parameters := + +;post-filter type := Gamma +;Gamma Filter Parameters := +;End Gamma Filter Parameters := + +output filename prefix := my_test_sim_image_SRT2DSPECT + + +end:= diff --git a/recon_test_pack/SRT2D_test_sim.par b/recon_test_pack/SRT2D_test_sim.par new file mode 100644 index 0000000000..fd5fed1cf8 --- /dev/null +++ b/recon_test_pack/SRT2D_test_sim.par @@ -0,0 +1,19 @@ +SRT2DParameters := +; test file for SRT2D +input file := my_precorrected_sino${suffix}.hs + +xy output image size (in pixels) := 91 +zoom := .5 + +;post-filter type := Wiener +;Wiener Filter Parameters := +;End Wiener Filter Parameters := + +;post-filter type := Gamma +;Gamma Filter Parameters := +;End Gamma Filter Parameters := + +output filename prefix := my_test_sim_image_SRT2D + + +end:= diff --git a/recon_test_pack/forward_project_SPECTUB.par b/recon_test_pack/forward_project_SPECTUB.par new file mode 100644 index 0000000000..ddcd2d6d6c --- /dev/null +++ b/recon_test_pack/forward_project_SPECTUB.par @@ -0,0 +1,35 @@ +Forward Projector parameters:= +; example par file for specifying the forward projector for e.g. fwdtest +type:=Matrix + Forward Projector Using Matrix Parameters := + + Matrix type := SPECT UB + Projection Matrix By Bin SPECT UB Parameters:= + ; same parameters as used in the OSMAPOSL_osem_SPECT.par file + + psf type:= Geometrical + + ;psf type:= 2D + ;maximum number of sigmas:= 2.0 + ;; sigma_at_depth = collimator_slope * depth_in_cm + collimator sigma 0(cm) + ;collimator slope := 0.0163 + ;collimator sigma 0(cm) := 0.1466 + + ;Attenuation correction { Simple // Full // No } + ;attenuation type := No + ;attenuation type := Simple + ;Values in attenuation map in cm-1 (float file) + ;attenuation map := attenuation.hv + attenuation map := my_atten_image_SPECT_modified.hv + + ;Mask properties { Cylinder // Attenuation Map // Explicit Mask // No} + ;mask type := Attenuation Map + mask type := No + + End Projection Matrix By Bin SPECT UB Parameters:= + + End Forward Projector Using Matrix Parameters := +end:= + + + diff --git a/recon_test_pack/generate_atten_cylinder_SPECT.par b/recon_test_pack/generate_atten_cylinder_SPECT.par new file mode 100644 index 0000000000..2f88210255 --- /dev/null +++ b/recon_test_pack/generate_atten_cylinder_SPECT.par @@ -0,0 +1,24 @@ +generate_image Parameters := +output filename:=my_atten_image_SPECT +X output image size (in pixels):=111 +Y output image size (in pixels):=111 +Z output image size (in pixels):=47 +X voxel size (in mm):= 3 +Y voxel size (in mm):= 3 +Z voxel size (in mm) :=3.27 + + Z number of samples to take per voxel := 1 + Y number of samples to take per voxel := 1 + X number of samples to take per voxel := 1 + +shape type:= ellipsoidal cylinder +Ellipsoidal Cylinder Parameters:= + radius-x (in mm):=100 + radius-y (in mm):=100 + length-z (in mm):=110 + origin (in mm):={70,10,-20} +;{70,10,-20} + END:= +value := 0.096 + +END:= diff --git a/recon_test_pack/run_scatter_tests.sh b/recon_test_pack/run_scatter_tests.sh index 8dd29343bb..237e52c339 100755 --- a/recon_test_pack/run_scatter_tests.sh +++ b/recon_test_pack/run_scatter_tests.sh @@ -37,7 +37,7 @@ examples_dir=`stir_config --examples-dir` scatter_pardir="$examples_dir/samples/scatter_estimation_par_files" echo "Using scatter parameter files from $scatter_pardir" -./simulate_PET_data_for_tests.sh +./simulate_data_for_tests.sh if [ $? -ne 0 ]; then echo "Error running simulation" exit 1 diff --git a/recon_test_pack/run_test_simulate_and_recon.sh b/recon_test_pack/run_test_simulate_and_recon.sh index 22c578a263..d5a282097d 100755 --- a/recon_test_pack/run_test_simulate_and_recon.sh +++ b/recon_test_pack/run_test_simulate_and_recon.sh @@ -9,9 +9,9 @@ # SPDX-License-Identifier: Apache-2.0 # # See STIR/LICENSE.txt for details -# +# # Author Kris Thielemans -# +# Author Dimitra Kyriakopoulou echo This script should work with STIR version 6.2. If you have echo a later version, you might have to update your test pack. @@ -58,27 +58,37 @@ fi echo "Using `command -v OSMAPOSL`" echo "Using `command -v OSSPS`" echo "Using `command -v FBP2D`" +echo "Using `command -v FBP3DRP`" +echo "Using `command -v SRT2D`" +echo "Using `command -v SRT2DSPECT`" # first need to set this to the C locale, as this is what the STIR utilities use # otherwise, awk might interpret floating point numbers incorrectly LC_ALL=C export LC_ALL -./simulate_PET_data_for_tests.sh +./simulate_data_for_tests.sh if [ $? -ne 0 ]; then echo "Error running simulation" exit 1 fi # need to repeat with zero-offset now as FBP doesn't support it zero_view_suffix=_force_zero_view_offset -./simulate_PET_data_for_tests.sh --force_zero_view_offset --suffix $zero_view_suffix +./simulate_data_for_tests.sh --force_zero_view_offset --suffix $zero_view_suffix if [ $? -ne 0 ]; then echo "Error running simulation with zero view offset" exit 1 fi ## TOF data TOF_suffix=_TOF -./simulate_PET_data_for_tests.sh --TOF --suffix "$TOF_suffix" +./simulate_data_for_tests.sh --TOF --suffix "$TOF_suffix" +if [ $? -ne 0 ]; then + echo "Error running simulation" + exit 1 +fi +## SPECT data +SPECT_suffix=_SPECT +./simulate_data_for_tests.sh --SPECT --suffix "$SPECT_suffix" if [ $? -ne 0 ]; then echo "Error running simulation" exit 1 @@ -96,7 +106,7 @@ input_ROI_mean=`awk 'NR>2 {print $2}' ${input_image}.roistats` # warning: currently OSMAPOSL needs to be run before OSSPS as # the OSSPS par file uses an OSMAPOSL result as initial image # and reuses its subset sensitivities -for recon in FBP2D FBP3DRP OSMAPOSL OSSPS; do +for recon in FBP2D FBP3DRP SRT2D SRT2DSPECT OSMAPOSL OSSPS ; do echo "========== Testing `command -v ${recon}`" # Check if we have CUDA code and parallelproj. # If so, check for test files in CUDA/* @@ -118,22 +128,31 @@ for recon in FBP2D FBP3DRP OSMAPOSL OSSPS; do for dataSuffix in "" "$TOF_suffix"; do echo "===== data suffix: \"$dataSuffix\"" # test first if analytic reconstruction and if so, run pre-correction - isFBP=0 + is_analytic=0 if expr "$recon" : FBP > /dev/null; then - if expr "$dataSuffix" : '.*TOF.*' > /dev/null; then - echo "Skipping TOF as not yet supported for FBP" - break - fi - isFBP=1 - suffix=$zero_view_suffix - export suffix - echo "Running precorrection" - correct_projdata correct_projdata_simulation.par > my_correct_projdata_simulation.log 2>&1 - if [ $? -ne 0 ]; then - echo "Error running precorrection. CHECK my_correct_projdata_simulation.log" - error_log_files="${error_log_files} my_correct_projdata_simulation.log" + is_analytic=1 + elif expr "$recon" : SRT > /dev/null; then + is_analytic=1 + fi + if [ $is_analytic = 1 ]; then + if expr "$dataSuffix" : '.*TOF.*' > /dev/null; then + echo "Skipping TOF as not yet supported for FBP and SRT" break - fi + fi + if expr "$recon" : SRT2DSPECT > /dev/null; then + suffix=$SPECT_suffix + export suffix + else + suffix=$zero_view_suffix + export suffix + echo "Running precorrection" + correct_projdata correct_projdata_simulation.par > my_correct_projdata_simulation.log 2>&1 + if [ $? -ne 0 ]; then + echo "Error running precorrection. CHECK my_correct_projdata_simulation.log" + error_log_files="${error_log_files} my_correct_projdata_simulation.log" + break + fi + fi else suffix="$dataSuffix" export suffix @@ -160,7 +179,7 @@ for recon in FBP2D FBP3DRP OSMAPOSL OSSPS; do output_filename=`awk -F':=' '/output[ _]*filename[ _]*prefix/ { value=$2;gsub(/[ \t]/, "", value); printf("%s", value) }' "$parfile"` # substitute env variables (e.g. to fill in suffix) output_filename=`eval echo "${output_filename}"` - if [ ${isFBP} -eq 0 ]; then + if [ ${is_analytic} -eq 0 ]; then # iterative algorithm, so we need to append the num_subiterations num_subiterations=`awk -F':=' '/number[ _]*of[ _]*subiterations/ { value=$2;gsub(/[ \t]/, "", value); printf("%s", value) }' ${parfile}` output_filename=${output_filename}_${num_subiterations} @@ -202,4 +221,3 @@ else tail ${error_log_files} exit 1 fi - diff --git a/recon_test_pack/simulate_data.sh b/recon_test_pack/simulate_data.sh index d8ed162eba..440d43cf80 100755 --- a/recon_test_pack/simulate_data.sh +++ b/recon_test_pack/simulate_data.sh @@ -3,6 +3,7 @@ # # Copyright (C) 2011 - 2011-01-14, Hammersmith Imanet Ltd # Copyright (C) 2011-07-01 - 2011, Kris Thielemans +# Copyright (C) 2024 University College London # This file is part of STIR. # # SPDX-License-Identifier: Apache-2.0 @@ -10,7 +11,10 @@ # See STIR/LICENSE.txt for details # # Author Kris Thielemans -# +# Author Dimitra Kyriakopoulou + +# Set single-threaded execution +#export OMP_NUM_THREADS=1 if [ $# -lt 3 -o $# -gt 5 ]; then echo "Usage: `basename $0` emission_image attenuation_image template_sino [ background_value [suffix] ]" @@ -29,51 +33,71 @@ fi if [ $# -gt 4 ]; then suffix=$5 fi -echo "=== create ACFs" -calculate_attenuation_coefficients --ACF my_acfs$suffix.hs ${atten_image} ${template_sino} forward_projector_proj_matrix_ray_tracing.par > my_create_acfs${suffix}.log 2>&1 -if [ $? -ne 0 ]; then - echo "ERROR running calculate_attenuation_coefficients. Check my_create_acfs${suffix}.log"; exit 1; -fi -echo "=== create line integrals" -forward_project my_line_integrals$suffix.hs ${emission_image} ${template_sino} forward_projector_proj_matrix_ray_tracing.par > my_create_line_integrals${suffix}.log 2>&1 -if [ $? -ne 0 ]; then - echo "ERROR running forward_project. Check my_create_line_integrals${suffix}.log"; exit 1; -fi +# Check if the suffix is "_SPECT" +if [ "$suffix" = "_SPECT" ]; then -echo "=== create constant randoms background" -${INSTALL_DIR}stir_math -s --including-first \ - --times-scalar 0 --add-scalar $background_value \ - my_randoms$suffix my_line_integrals$suffix.hs -if [ $? -ne 0 ]; then - echo "ERROR running stir_math"; exit 1; -fi + echo "=== create line integrals for SPECT sinogram" + forward_project "my_sino${suffix}.hs" "${emission_image}" "${template_sino}" forward_project_SPECTUB.par > "my_create_sino${suffix}.log" 2>&1 + if [ $? -ne 0 ]; then + echo "ERROR running forward_project for SPECT sinogram. Check my_create_sino${suffix}.log"; exit 1; + fi -echo "=== create norm factors" -# currently just 1 as not used in rest of script yet. -stir_math -s --including-first \ - --times-scalar 0 --add-scalar 1 my_norm$suffix.hs my_acfs$suffix.hs - -echo "=== create prompts" -INPUT=my_line_integrals${suffix}.hs OUTPUT=my_prompts${suffix}.hs \ - MULT=my_acfs${suffix}.hs \ - RANDOMS=my_randoms${suffix}.hs \ - correct_projdata uncorrect_projdata.par > my_create_prompts${suffix}.log 2>&1 -if [ $? -ne 0 ]; then - echo "ERROR running correct_projdata. Check my_create_prompts${suffix}.log"; exit 1; -fi -# could call poisson_noise here + echo "=== create line integrals for SRT2DSPECT attenuation projection" + forward_project "my_attenuation_sino$suffix.hs" "${atten_image}" "${template_sino}" > "my_create_attenuation_sino${suffix}.log" 2>&1 + if [ $? -ne 0 ]; then + echo "ERROR running forward_project for attenuation sinogram. Check my_create_attenuation_sino${suffix}.log"; exit 1; + fi + +else + + echo "=== create ACFs" + calculate_attenuation_coefficients --ACF my_acfs$suffix.hs ${atten_image} ${template_sino} forward_projector_proj_matrix_ray_tracing.par > my_create_acfs${suffix}.log 2>&1 + if [ $? -ne 0 ]; then + echo "ERROR running calculate_attenuation_coefficients. Check my_create_acfs${suffix}.log"; exit 1; + fi + + echo "=== create line integrals" + forward_project my_line_integrals$suffix.hs ${emission_image} ${template_sino} forward_projector_proj_matrix_ray_tracing.par > my_create_line_integrals${suffix}.log 2>&1 + if [ $? -ne 0 ]; then + echo "ERROR running forward_project. Check my_create_line_integrals${suffix}.log"; exit 1; + fi + + echo "=== create constant randoms background" + ${INSTALL_DIR}stir_math -s --including-first \ + --times-scalar 0 --add-scalar $background_value \ + my_randoms$suffix my_line_integrals$suffix.hs + if [ $? -ne 0 ]; then + echo "ERROR running stir_math"; exit 1; + fi + + echo "=== create norm factors" + # currently just 1 as not used in rest of script yet. + stir_math -s --including-first \ + --times-scalar 0 --add-scalar 1 my_norm$suffix.hs my_acfs$suffix.hs + + echo "=== create prompts" + INPUT=my_line_integrals${suffix}.hs OUTPUT=my_prompts${suffix}.hs \ + MULT=my_acfs${suffix}.hs \ + RANDOMS=my_randoms${suffix}.hs \ + correct_projdata uncorrect_projdata.par > my_create_prompts${suffix}.log 2>&1 + if [ $? -ne 0 ]; then + echo "ERROR running correct_projdata. Check my_create_prompts${suffix}.log"; exit 1; + fi + + # could call poisson_noise here -echo "=== create additive sinogram for reconstruction" -# need randoms (and scatter) multiplied by ACF and norm (but we don't have a norm here yet) -# need to use correct_projdata as in TOF, ACF/norm is non-TOF, so stir_math will fail -INPUT=my_randoms${suffix}.hs OUTPUT=my_additive_sinogram${suffix}.hs \ -MULT=my_acfs${suffix}.hs \ - correct_projdata correct_projdata_norm_only.par > my_create_additive_sino${suffix}.log 2>&1 -if [ $? -ne 0 ]; then - echo "ERROR running correct_projdata. Check my_create_additive_sino${suffix}.log"; exit 1; + echo "=== create additive sinogram for reconstruction" + # need randoms (and scatter) multiplied by ACF and norm (but we don't have a norm here yet) + # need to use correct_projdata as in TOF, ACF/norm is non-TOF, so stir_math will fail + INPUT=my_randoms${suffix}.hs OUTPUT=my_additive_sinogram${suffix}.hs \ + MULT=my_acfs${suffix}.hs \ + correct_projdata correct_projdata_norm_only.par > my_create_additive_sino${suffix}.log 2>&1 + if [ $? -ne 0 ]; then + echo "ERROR running correct_projdata. Check my_create_additive_sino${suffix}.log"; exit 1; + fi fi echo "Done creating simulated data" diff --git a/recon_test_pack/simulate_PET_data_for_tests.sh b/recon_test_pack/simulate_data_for_tests.sh similarity index 61% rename from recon_test_pack/simulate_PET_data_for_tests.sh rename to recon_test_pack/simulate_data_for_tests.sh index 6f31b26189..18e04f42eb 100755 --- a/recon_test_pack/simulate_PET_data_for_tests.sh +++ b/recon_test_pack/simulate_data_for_tests.sh @@ -15,7 +15,7 @@ # See STIR/LICENSE.txt for details # # Author Kris Thielemans -# +# Author Dimitra Kyriakopoulou echo This script should work with STIR version 6.x. If you have echo a later version, you might have to update your test pack. @@ -29,6 +29,7 @@ generate_images=1 force_zero_view_offset=0 TOF=0 suffix="" +SPECT=0 # # Parse option arguments (--) # Note that the -- is required to suppress interpretation of $1 as options @@ -43,6 +44,9 @@ do elif test "$1" = "--TOF" then TOF=1 + elif test "$1" = "--SPECT" + then + SPECT=1 elif test "$1" = "--suffix" then suffix="$2" @@ -61,7 +65,7 @@ do echo "Randoms+scatter are set to a constant proj_data. you can set that value via the environment variable:" echo " background_value" exit 1 - else +else echo Warning: Unknown option "$1" echo rerun with --help for more info. exit 1 @@ -86,8 +90,34 @@ if [ "$generate_images" -eq 1 ]; then generate_image generate_uniform_cylinder.par echo "=== make attenuation image" generate_image generate_atten_cylinder.par +echo "=== make attenuation image for SPECT: to make up for the fliprl of SPECTUB" +generate_image generate_atten_cylinder_SPECT.par fi +# Function to comment out the specific line in a given file and write to a new file +comment_out_line() { + input_file="$1" + output_file="$2" + + sed '2s/^/;/' "$input_file" > "$output_file" +} + +# Paths to the input and output files for uniform cylinder +uniform_input_file="my_uniform_cylinder.hv" +uniform_output_file="my_uniform_cylinder_SPECT.hv" + +# Comment out the "modality" line in the uniform image file +comment_out_line "$uniform_input_file" "$uniform_output_file" + +# Paths to the input and output files for attenuation image +atten_input_file="my_atten_image_SPECT.hv" +atten_output_file="my_atten_image_SPECT_modified.hv" + +# Comment out the specific line in the attenuation image file +comment_out_line "$atten_input_file" "$atten_output_file" + + +if [ "$SPECT" -eq 0 ]; then if [ "$TOF" -eq 0 ]; then : ${view_mash:=1} : ${span:=2} @@ -121,31 +151,42 @@ $span $max_rd EOF fi + create_projdata_template ${template_sino} < my_input.txt > my_create_${template_sino}.log 2>&1 if [ $? -ne 0 ]; then echo "ERROR running create_projdata_template. Check my_create_${template_sino}.log"; exit 1; fi -# fix-up header by insert energy info just before the end -# trick for awk comes from the www.theunixschool.com -awk '/END OF INTERFILE/ { print "number of energy windows := 1\nenergy window lower level[1] := 350\nenergy window upper level[1] := 650\nEnergy resolution := 0.22\nReference energy (in keV) := 511" }1 ' \ - ${template_sino} > tmp_header.hs -mv tmp_header.hs ${template_sino} - -if [ $force_zero_view_offset -eq 1 ]; then - if [ "$TOF" -eq 1 ]; then - echo "$0 would need work to be used with both TOF and zero-offset. Exiting" - exit 1 + # fix-up header by insert energy info just before the end + # trick for awk comes from the www.theunixschool.com + awk '/END OF INTERFILE/ { print "number of energy windows := 1\nenergy window lower level[1] := 350\nenergy window upper level[1] := 650\nEnergy resolution := 0.22\nReference energy (in keV) := 511" }1 ' \ + "${template_sino}" > tmp_header.hs + mv tmp_header.hs "${template_sino}" + if [ $force_zero_view_offset -eq 1 ]; then + if [ "$TOF" -eq 1 ]; then + echo "$0 would need work to be used with both TOF and zero-offset. Exiting" + exit 1 + fi + new_template_sino="my_DSTE_3D_rd2_template${suffix}.hs" + force_view_offset_to_zero.sh "${new_template_sino}" "${template_sino}" + template_sino="${new_template_sino}" fi - new_template_sino=my_DSTE_3D_rd2_template$suffix.hs - force_view_offset_to_zero.sh ${new_template_sino} ${template_sino} - template_sino=${new_template_sino} fi + # create sinograms : ${background_value:=10} -./simulate_data.sh my_uniform_cylinder.hv my_atten_image.hv ${template_sino} ${background_value} ${suffix} -if [ $? -ne 0 ]; then - echo "Error running simulation" - exit 1 +if [ "$SPECT" -eq 1 ]; then + # create SPECT sinograms + ./simulate_data.sh "my_uniform_cylinder_SPECT.hv" "my_atten_image_SPECT_modified.hv" "SPECT_test_Interfile_header.hs" "${background_value}" "${suffix}" + if [ $? -ne 0 ]; then + echo "Error running SPECT simulation" + exit 1 + fi +else + ./simulate_data.sh "my_uniform_cylinder.hv" "my_atten_image.hv" "${template_sino}" "${background_value}" "${suffix}" + if [ $? -ne 0 ]; then + echo "Error running PET simulation" + exit 1 + fi fi diff --git a/scripts/maintenance/make_distribution.sh b/scripts/maintenance/make_distribution.sh index 0f95f9219e..10eeacd39a 100755 --- a/scripts/maintenance/make_distribution.sh +++ b/scripts/maintenance/make_distribution.sh @@ -64,7 +64,7 @@ set -e #fi -read -p "Did you update CMakeLists.txt, version numbers in \*tex files, documentation/history.htm, CITATION.cff? (press Ctrl-C if not)" +read -p "Did you update CMakeLists.txt, version numbers in \*tex files, documentation/history.htm, documentation/credits.htm, CITATION.cff? (press Ctrl-C if not)" mkdir -p ${DISTRIB} cd ${DISTRIB} @@ -257,4 +257,4 @@ if [ $do_website_final_version = 0 ]; then echo "if not beta, did you run with 'do_website_final_version=1'?" fi echo "create a GitHub release. To get release notes, use" -echo " pandoc -t markdown_github -o release_${VERSION%%.?}.md release_${VERSION%%.?}.htm" +echo " pandoc -t gfm -o release_${VERSION%%.?}.md STIR/documentation/release_${VERSION%%.?}.htm" diff --git a/src/IO/CMakeLists.txt b/src/IO/CMakeLists.txt index ef8e53e6da..90f4115f40 100644 --- a/src/IO/CMakeLists.txt +++ b/src/IO/CMakeLists.txt @@ -83,9 +83,9 @@ endif() if (CERN_ROOT_FOUND) target_include_directories(IO PRIVATE ${CERN_ROOT_INCLUDE_DIRS}) if (TARGET ROOT::Tree) - target_link_libraries(IO PUBLIC ROOT::Tree) + target_link_libraries(IO PRIVATE ROOT::Tree) else() - target_link_libraries(IO PUBLIC ${CERN_ROOT_LIBRARIES}) + target_link_libraries(IO PRIVATE ${CERN_ROOT_LIBRARIES}) endif() endif() @@ -94,12 +94,12 @@ if (HAVE_HDF5) endif() if (HAVE_ITK) - target_link_libraries(IO PUBLIC ITKCommon ${ITK_LIBRARIES}) + target_link_libraries(IO PRIVATE ITKCommon ${ITK_LIBRARIES}) endif() if (UPENN_FOUND) target_include_directories(IO PUBLIC ${UPENN_INCLUDE_DIR}) - target_link_libraries(IO PUBLIC ${UPENN_libsss_tof} ${UPENN_libfit} ${UPENN_libdist} ${UPENN_libgeom} + target_link_libraries(IO PRIVATE ${UPENN_libsss_tof} ${UPENN_libfit} ${UPENN_libdist} ${UPENN_libgeom} ${UPENN_liblor} ${UPENN_liblist} ${UPENN_libmhdr} ${JANSSON_LIBRARY} ${ZLIB_LIBRARY_RELEASE} ${UPENN_libimagio} ${UPENN_libimagio++}) endif() diff --git a/src/IO/GEHDF5Wrapper.cxx b/src/IO/GEHDF5Wrapper.cxx index 45978ba2ed..bc6abad67c 100644 --- a/src/IO/GEHDF5Wrapper.cxx +++ b/src/IO/GEHDF5Wrapper.cxx @@ -65,6 +65,17 @@ read_float(const H5::H5File& file, const std::string& dataset) return tmp; } +static std::string +read_string(const H5::H5File& file, const std::string& dataset) +{ + H5::DataSet ds = file.openDataSet(dataset.c_str()); + H5::StrType datatype = ds.getStrType(); + + std::string value; + ds.read(value, datatype); + return value; +} + bool GEHDF5Wrapper::check_GE_signature(const std::string& filename) { @@ -90,13 +101,7 @@ GEHDF5Wrapper::check_GE_signature(H5::H5File& file) if (file.getId() == -1) error("File is not open. Aborting"); - H5::StrType vlst( - 0, - 37); // 37 here is the length of the string (PW got it from the text file generated by list2txt with the LIST000_decomp.BLF - std::string read_str_manufacturer; - - H5::DataSet dataset2 = file.openDataSet("/HeaderData/ExamData/manufacturer"); - dataset2.read(read_str_manufacturer, vlst); + std::string read_str_manufacturer = read_string(file, "/HeaderData/ExamData/manufacturer"); if (read_str_manufacturer == "GE MEDICAL SYSTEMS") { @@ -277,14 +282,7 @@ GEHDF5Wrapper::open(const std::string& filename) shared_ptr GEHDF5Wrapper::get_scanner_from_HDF5() { - std::string read_str_scanner; - H5::StrType vlst( - 0, - 37); // 37 here is the length of the string (PW got it from the text file generated by list2txt with the LIST000_decomp.BLF - - H5::DataSet dataset = file.openDataSet("/HeaderData/ExamData/scannerDesc"); - dataset.read(read_str_scanner, vlst); - + std::string read_str_scanner = read_string(this->file, "/HeaderData/ExamData/scannerDesc"); float effective_ring_diameter; int num_transaxial_blocks_per_bucket = 0; int num_axial_blocks_per_bucket = 0; @@ -425,6 +423,11 @@ GEHDF5Wrapper::initialise_proj_data_info_from_HDF5() { shared_ptr scanner_sptr = get_scanner_from_HDF5(); + // TODO get TOF mashing when reading sinos as TOF + const auto num_tof_bins = read_dataset_uint32("/HeaderData/Sorter/numTOF_bins"); + if (num_tof_bins > 1) + warning("GE RDF data currently still read as non-TOF"); + this->proj_data_info_sptr = ProjDataInfo::construct_proj_data_info(scanner_sptr, /*span*/ 2, @@ -654,11 +657,20 @@ GEHDF5Wrapper::initialise_proj_data(const unsigned int view_num) if (view_num == 0 || view_num > static_cast(this->get_scanner_sptr()->get_num_detectors_per_ring() / 2)) error("internal error in GE HDF5 code: view number " + std::to_string(view_num) + " is incorrect"); + const auto num_tof_bins = read_dataset_uint32("/HeaderData/Sorter/numTOF_bins"); + if (rdf_ver == 9) { - m_address = "/SegmentData/Segment2/3D_TOF_Sinogram/view" + std::to_string(view_num); - - m_dataset_sptr.reset(new H5::DataSet(file.openDataSet(m_address))); + if (num_tof_bins > 1) + { + m_address = "/SegmentData/Segment2/3D_TOF_Sinogram/view" + std::to_string(view_num); + m_dataset_sptr.reset(new H5::DataSet(file.openDataSet(m_address))); + } + else + { + m_address = "/SegmentData/Segment2/3D_Sinogram/view" + std::to_string(view_num); + m_dataset_sptr.reset(new H5::DataSet(file.openDataSet(m_address))); + } m_dataspace = m_dataset_sptr->getSpace(); // Create an array to host the size of the dimensions const int rank = m_dataspace.getSimpleExtentNdims(); @@ -673,13 +685,7 @@ GEHDF5Wrapper::initialise_proj_data(const unsigned int view_num) // AB for signa, these where [1981,27,357] and [45,448,357] m_NX_SUB = dims[0]; // hyperslab dimensions m_NY_SUB = dims[1]; - m_NZ_SUB = dims[2]; -#if 0 - // AB todo: ??? why are these different? - m_NX = 45; // output buffer dimensions - m_NY = 448; - m_NZ = 357; -#endif + m_NZ_SUB = rank > 2 ? dims[2] : 1; } else return Succeeded::no; @@ -831,7 +837,6 @@ GEHDF5Wrapper::read_sinogram(Array<3, unsigned char>& output, } } } - output.release_data_ptr(); } return Succeeded::yes; diff --git a/src/IO/MultiDynamicDiscretisedDensityOutputFileFormat.cxx b/src/IO/MultiDynamicDiscretisedDensityOutputFileFormat.cxx index 4542fcabd5..7ee8c09dc1 100644 --- a/src/IO/MultiDynamicDiscretisedDensityOutputFileFormat.cxx +++ b/src/IO/MultiDynamicDiscretisedDensityOutputFileFormat.cxx @@ -93,7 +93,7 @@ MultiDynamicDiscretisedDensityOutputFileFormat::actual_write_to_file(std::string // Create all the filenames VectorWithOffset individual_filenames(1, int(density.get_num_time_frames())); for (int i = 1; i <= int(density.get_num_time_frames()); i++) - individual_filenames[i] = filename + "_" + boost::lexical_cast(i); + individual_filenames[i] = filename + "_" + std::to_string(i); // Write each individual image for (int i = 1; i <= int(density.get_num_time_frames()); i++) diff --git a/src/IO/MultiParametricDiscretisedDensityOutputFileFormat.cxx b/src/IO/MultiParametricDiscretisedDensityOutputFileFormat.cxx index 126c5bd5f0..b065304d85 100644 --- a/src/IO/MultiParametricDiscretisedDensityOutputFileFormat.cxx +++ b/src/IO/MultiParametricDiscretisedDensityOutputFileFormat.cxx @@ -106,7 +106,7 @@ ParamDiscDensityOutputFileFormat::actual_write_to_file(std::string& filename, // Create all the filenames VectorWithOffset individual_filenames(1, int(density.get_num_params())); for (int i = 1; i <= int(density.get_num_params()); i++) - individual_filenames[i] = filename + "_" + boost::lexical_cast(i); + individual_filenames[i] = filename + "_" + std::to_string(i); // Write each individual image for (int i = 1; i <= int(density.get_num_params()); i++) diff --git a/src/analytic/SRT2D/CMakeLists.txt b/src/analytic/SRT2D/CMakeLists.txt new file mode 100644 index 0000000000..587031b5c8 --- /dev/null +++ b/src/analytic/SRT2D/CMakeLists.txt @@ -0,0 +1,23 @@ +# +# + +set(dir analytic_SRT2D) + +set (dir_LIB_SOURCES ${dir}_LIB_SOURCES) + +set(${dir_LIB_SOURCES} + SRT2DReconstruction.cxx +) + +#$(dir)_REGISTRY_SOURCES:= + +include(stir_lib_target) +target_link_libraries(analytic_SRT2D recon_buildblock IO ) + +set (dir_EXE_SOURCES ${dir}_EXE_SOURCES) + +set(${dir_EXE_SOURCES} + SRT2D.cxx +) + +include(stir_exe_targets) diff --git a/src/analytic/SRT2D/SRT2D.cxx b/src/analytic/SRT2D/SRT2D.cxx new file mode 100644 index 0000000000..844b0aa106 --- /dev/null +++ b/src/analytic/SRT2D/SRT2D.cxx @@ -0,0 +1,36 @@ +// +// +/*! + \file + \ingroup analytic + \brief Main program for SRT2D reconstruction + \author Dimitra Kyriakopoulou + \author Kris Thielemans +*/ +/* + Copyright (C) 2024, University College London + + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +#include "stir/analytic/SRT2D/SRT2DReconstruction.h" +#include "stir/Succeeded.h" +#ifndef PARALLEL +# define Main main +#else +# define Main master_main +#endif + +USING_NAMESPACE_STIR + +int +Main(int argc, char** argv) +{ + SRT2DReconstruction reconstruction_object(argc > 1 ? argv[1] : ""); + + return reconstruction_object.reconstruct() == Succeeded::yes ? EXIT_SUCCESS : EXIT_FAILURE; +} diff --git a/src/analytic/SRT2D/SRT2DReconstruction.cxx b/src/analytic/SRT2D/SRT2DReconstruction.cxx new file mode 100644 index 0000000000..60ffbd14a4 --- /dev/null +++ b/src/analytic/SRT2D/SRT2DReconstruction.cxx @@ -0,0 +1,604 @@ +/* + Copyright (C) 2024 University College London + + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ +/*! + \file + \ingroup analytic + \brief Implementation of class stir::SRT2DReconstruction + + \author Dimitra Kyriakopoulou + \author Kris Thielemans +*/ + +#include "stir/analytic/SRT2D/SRT2DReconstruction.h" +#include "stir/VoxelsOnCartesianGrid.h" +#include "stir/RelatedViewgrams.h" +#include "stir/recon_buildblock/BackProjectorByBinUsingInterpolation.h" +#include "stir/ProjDataInfoCylindricalArcCorr.h" +#include "stir/ArcCorrection.h" +#include "stir/SSRB.h" +#include "stir/ProjDataInMemory.h" +#include "stir/Bin.h" +#include "stir/round.h" +#include "stir/display.h" +#include +#include "stir/IO/interfile.h" +#include "stir/info.h" +#include + +#ifdef STIR_OPENMP +# include +#endif + +#include // For M_PI and other math functions +#ifndef M_PI +# define M_PI 3.14159265358979323846 +#endif + +#ifdef STIR_OPENMP +# include "stir/num_threads.h" +#endif + +#include +#include + +START_NAMESPACE_STIR + +const char* const SRT2DReconstruction::registered_name = "SRT2D"; + +void +SRT2DReconstruction::set_defaults() +{ + base_type::set_defaults(); + num_segments_to_combine = -1; +} + +void +SRT2DReconstruction::initialise_keymap() +{ + base_type::initialise_keymap(); + + parser.add_start_key("SRT2DParameters"); + parser.add_stop_key("End"); + parser.add_key("num_segments_to_combine with SSRB", &num_segments_to_combine); +} + +void +SRT2DReconstruction::ask_parameters() +{ + base_type::ask_parameters(); + num_segments_to_combine = ask_num("num_segments_to_combine (must be odd)", -1, 101, -1); +} + +bool +SRT2DReconstruction::post_processing() +{ + return base_type::post_processing(); +} + +Succeeded +SRT2DReconstruction::set_up(shared_ptr const& target_data_sptr) +{ + if (base_type::set_up(target_data_sptr) == Succeeded::no) + return Succeeded::no; + + if (num_segments_to_combine >= 0 && num_segments_to_combine % 2 == 0) + error(boost::format("num_segments_to_combine has to be odd (or -1), but is %d") % num_segments_to_combine); + + if (num_segments_to_combine == -1) + { + const shared_ptr proj_data_info_cyl_sptr + = dynamic_pointer_cast(proj_data_ptr->get_proj_data_info_sptr()); + + if (is_null_ptr(proj_data_info_cyl_sptr)) + num_segments_to_combine = 1; // cannot SSRB non-cylindrical data yet + else + { + if (proj_data_info_cyl_sptr->get_min_ring_difference(0) != proj_data_info_cyl_sptr->get_max_ring_difference(0) + || proj_data_info_cyl_sptr->get_num_segments() == 1) + num_segments_to_combine = 1; + else + num_segments_to_combine = 3; + } + } + return Succeeded::yes; +} + +std::string +SRT2DReconstruction::method_info() const +{ + return "SRT2D"; +} + +SRT2DReconstruction::SRT2DReconstruction(const std::string& parameter_filename) +{ + initialise(parameter_filename); + info(boost::format("%1%") % parameter_info()); +} + +SRT2DReconstruction::SRT2DReconstruction() +{ + set_defaults(); +} + +SRT2DReconstruction::SRT2DReconstruction(const shared_ptr& proj_data_ptr_v, const int num_segments_to_combine_v) +{ + set_defaults(); + proj_data_ptr = proj_data_ptr_v; + num_segments_to_combine = num_segments_to_combine_v; +} + +Succeeded +SRT2DReconstruction::actual_reconstruct(shared_ptr> const& density_ptr) +{ + // In case of 3D data, use only direct sinograms + // perform SSRB + if (num_segments_to_combine > 1) + { + info(boost::format("Performing SSRB with num_segments_to_combine = %d") % num_segments_to_combine); + const ProjDataInfoCylindrical& proj_data_info_cyl + = dynamic_cast(*proj_data_ptr->get_proj_data_info_sptr()); + + shared_ptr ssrb_info_sptr( + SSRB(proj_data_info_cyl, num_segments_to_combine, 1, 0, (num_segments_to_combine - 1) / 2)); + shared_ptr proj_data_to_SRT_ptr(new ProjDataInMemory(proj_data_ptr->get_exam_info_sptr(), ssrb_info_sptr)); + SSRB(*proj_data_to_SRT_ptr, *proj_data_ptr); + proj_data_ptr = proj_data_to_SRT_ptr; + } + else + { + info("Using existing proj_data_ptr"); + } + + // check if segment 0 has direct sinograms + { + const float tan_theta = proj_data_ptr->get_proj_data_info_sptr()->get_tantheta(Bin(0, 0, 0, 0)); + if (fabs(tan_theta) > 1.E-4) + { + warning("SRT2D: segment 0 has non-zero tan(theta) %g", tan_theta); + return Succeeded::no; + } + } + + float tangential_sampling; + shared_ptr arc_corrected_proj_data_info_sptr; + + // arc-correction if necessary + ArcCorrection arc_correction; + bool do_arc_correction = false; + if (!is_null_ptr(dynamic_pointer_cast(proj_data_ptr->get_proj_data_info_sptr()))) + { + info("Data is already arc-corrected"); + arc_corrected_proj_data_info_sptr = proj_data_ptr->get_proj_data_info_sptr()->create_shared_clone(); + tangential_sampling = dynamic_cast(*proj_data_ptr->get_proj_data_info_sptr()) + .get_tangential_sampling(); + } + else + { + info("Performing arc-correction"); + if (arc_correction.set_up(proj_data_ptr->get_proj_data_info_sptr()->create_shared_clone()) == Succeeded::no) + return Succeeded::no; + do_arc_correction = true; + arc_corrected_proj_data_info_sptr = arc_correction.get_arc_corrected_proj_data_info_sptr(); + tangential_sampling = arc_correction.get_arc_corrected_proj_data_info().get_tangential_sampling(); + } + + VoxelsOnCartesianGrid& image = dynamic_cast&>(*density_ptr); + Sinogram sino = proj_data_ptr->get_empty_sinogram(0, 0); + Viewgram view = proj_data_ptr->get_viewgram(0, 0); + if (do_arc_correction) + { + info("Arc-correcting viewgram"); + view = arc_correction.do_arc_correction(view); + } + Viewgram view1 = proj_data_ptr->get_empty_viewgram(0, 0); + Viewgram view_th = proj_data_ptr->get_empty_viewgram(0, 0); + Viewgram view1_th = proj_data_ptr->get_empty_viewgram(0, 0); + + // Retrieve runtime-dependent sizes + const int sp = view.get_num_tangential_poss(); + const int sth = proj_data_ptr->get_num_views(); + const int sa = proj_data_ptr->get_num_axial_poss(0); + const int sx = image.get_x_size(); + const int sy = image.get_y_size(); + const int sx2 = ceil(sx / 2.0), sy2 = ceil(sy / 2.0); + const int sth2 = ceil(sth / 2.0); + const float c_int = -1 / (M_PI * sth * (sp - 1)) * 2; // instead of *2 we will write /zoom; + + int ia, image_pos; + int ith, jth, ip, ix1, ix2, k1, k2; + + float x, aux; + + const int image_min_x = image.get_min_x(); + const int image_min_y = image.get_min_y(); + + // Dynamic declarations using std::vector + std::vector th(sth); // Theta values for each view angle. + std::vector p(sp); // Tangential positions for projections. + std::vector p_ud(sp); // Tangential positions for up-down flipping. + std::vector x1(sx); // X-coordinates for the reconstructed image grid. + std::vector x2(sy); // Y-coordinates for the reconstructed image grid. + + std::vector> f(sth, std::vector(sp, 0.0f)); // Projection data matrix. + std::vector> ddf(sth, std::vector(sp, 0.0f)); // Second derivative of projections. + + std::vector> f_ud(sth, std::vector(sp, 0.0f)); // Flipped projection data. + std::vector> ddf_ud(sth, std::vector(sp, 0.0f)); // Second derivatives of flipped projections. + + std::vector> f1(sth, std::vector(sp, 0.0f)); // Left-right mirrored projection data. + std::vector> ddf1(sth, + std::vector(sp, 0.0f)); // Second derivatives of left-right mirrored projections. + + std::vector> f1_ud(sth, std::vector(sp, 0.0f)); // Up-down flipped mirrored projection data. + std::vector> ddf1_ud( + sth, std::vector(sp, 0.0f)); // Second derivatives of up-down flipped mirrored projections. + + std::vector> f_th(sth, std::vector(sp, 0.0f)); // Projection data along theta views. + std::vector> ddf_th(sth, std::vector(sp, 0.0f)); // Second derivatives along theta views. + + std::vector> f_th_ud(sth, std::vector(sp, 0.0f)); // Flipped data along theta views. + std::vector> ddf_th_ud(sth, std::vector(sp, 0.0f)); // Second derivatives of flipped data. + + std::vector> f1_th(sth, std::vector(sp, 0.0f)); // Mirrored data along theta views. + std::vector> ddf1_th(sth, std::vector(sp, 0.0f)); // Second derivatives of mirrored data. + + std::vector> f1_th_ud(sth, std::vector(sp, 0.0f)); // Flipped mirrored data along theta views. + std::vector> ddf1_th_ud(sth, std::vector(sp, 0.0f)); // Second derivatives of flipped mirrored data. + + std::vector lg(sp); // Logarithmic differences for interpolation. + std::vector termC(sth); // Correction term for each view. + std::vector termC_th(sth); // Correction term for theta projections. + + const float dp6 = 6.0 / 4.0 * 2.0 / (sp - 1.0); // Integration constant for second derivatives. + +#ifdef STIR_OPENMP + set_num_threads(); +# pragma omp single + info("Using OpenMP-version of SRT2D with " + std::to_string(omp_get_num_threads()) + " threads on " + + std::to_string(omp_get_num_procs()) + " processors."); +#endif + + // Put theta and p in arrays. + for (ith = 0; ith < sth; ith++) + th[ith] = ith * M_PI / sth; + + for (ip = 0; ip < sp; ip++) + p[ip] = -1.0 + 2.0 * ip / (sp - 1); + for (ip = 0; ip < sp; ip++) + p_ud[sp - ip - 1] = p[ip]; + + //-- Creation of the grid + for (k1 = 0; k1 < sx; k1++) + x1[k1] = -1.0 * sx / (sp + 1) + 2.0 * sx / (sp + 1) * k1 / (sx - 1); + // x1[k1] = -1.0 * sx / ((sp + 1) * zoom) + k1 * 2.0 * sx / ((sp + 1) * zoom) / (sx - 1); + for (k2 = 0; k2 < sx; k2++) + x2[k2] = -1.0 * sx / (sp + 1) + 2.0 * sx / (sp + 1) * k2 / (sx - 1); + // x2[k2] = -1.0 * sx / ((sp + 1) * zoom) + k2 * 2.0 * sx / ((sp + 1) * zoom) / (sx - 1); + + // Starting calculations per view + // 2D algorithm only + + // ----- + // special case of ith=0 + // ----- + for (ia = 0; ia < sa; ia++) + { + for (ip = 0; ip < sp; ip++) + { + f[ia][ip] = view[view.get_min_axial_pos_num() + ia][view.get_min_tangential_pos_num() + ip]; + } + spline(p, f[ia], sp, ddf[ia]); + } + + for (ia = 0; ia < sa; ia++) + { + termC[ia] = (ddf[ia][0] * (3 * p[1] - p[0]) + ddf[ia][sp - 1] * (p[sp - 1] - 3.0 * p[sp - 2])) / 4.0; + for (ip = 0; ip < sp; ip++) + { + termC[ia] += dp6 * ddf[ia][ip]; + } + } + + for (ix1 = 0; ix1 < sx2; ix1++) + { + for (ix2 = 0; ix2 < sy2; ix2++) + { + aux = sqrt(1.0 - x2[ix2] * x2[ix2]); + if (fabs(x2[ix2]) >= 1.0 || fabs(x1[ix1]) >= aux) + { + continue; + } + x = x2[ix2] * cos(th[0]) - x1[ix1] * sin(th[0]); + + for (ip = 0; ip < sp; ip++) + { + double val = fabs(x - p[ip]); + lg[ip] = val < 2e-6 ? 0. : std::log(val); + } + + for (ia = 0; ia < sa; ia++) + { + int img_x = image_min_x + sx - ix1 - 1; + int img_y = image_min_y + ix2; + + image[ia][img_x][img_y] = -hilbert_der(x, f[ia], ddf[ia], p, sp, lg, termC[ia]) / (M_PI * sth * (sp - 1)); + } + } + } + +// jth, ia, ip, termC_th, ix2, ix1, aux, x, image_pos +#ifdef STIR_OPENMP +# pragma omp parallel firstprivate(f, \ + ddf, \ + f_ud, \ + ddf_ud, \ + f1, \ + ddf1, \ + f1_ud, \ + ddf1_ud, \ + f_th, \ + ddf_th, \ + f_th_ud, \ + ddf_th_ud, \ + f1_th, \ + ddf1_th, \ + f1_th_ud, \ + ddf1_th_ud, \ + termC, \ + termC_th, \ + lg) \ + shared(view, view1, view_th, view1_th, do_arc_correction, arc_correction, p_ud, p, th, x1, x2, image) private( \ + jth, ia, ip, ix2, ix1, aux, x, image_pos) +# pragma omp for schedule(static) nowait +#endif + for (ith = 1; ith < sth; ith++) + { + if (ith < sth2) + { + jth = sth2 - ith; + } + else if (ith > sth2) + { + jth = (int)ceil(3 * sth / 2.0) - ith; // MARK integer division + } + else + { + jth = sth2; + } + + // Loading related viewgrams +#ifdef STIR_OPENMP +# pragma omp critical +#endif + { + view = proj_data_ptr->get_viewgram(ith, 0); + view1 = proj_data_ptr->get_viewgram(sth - ith, 0); + view_th = proj_data_ptr->get_viewgram(jth, 0); + view1_th = proj_data_ptr->get_viewgram(sth - jth, 0); + if (do_arc_correction) + { + view = arc_correction.do_arc_correction(view); + view1 = arc_correction.do_arc_correction(view1); + view_th = arc_correction.do_arc_correction(view_th); + view1_th = arc_correction.do_arc_correction(view1_th); + } + + for (ia = 0; ia < sa; ia++) + { + for (ip = 0; ip < sp; ip++) + { + f[ia][ip] = view[view.get_min_axial_pos_num() + ia][view.get_min_tangential_pos_num() + ip]; + f_ud[ia][sp - ip - 1] = f[ia][ip]; + f1[ia][ip] = view1[view.get_min_axial_pos_num() + ia][view.get_min_tangential_pos_num() + ip]; + f1_ud[ia][sp - ip - 1] = f1[ia][ip]; + + f_th[ia][ip] = view_th[view.get_min_axial_pos_num() + ia][view.get_min_tangential_pos_num() + ip]; + f_th_ud[ia][sp - ip - 1] = f_th[ia][ip]; + f1_th[ia][ip] = view1_th[view.get_min_axial_pos_num() + ia][view.get_min_tangential_pos_num() + ip]; + f1_th_ud[ia][sp - ip - 1] = f1_th[ia][ip]; + } + } + } + + // Calculation of second derivative by use of function spline + for (ia = 0; ia < sa; ia++) + { + spline(p, f[ia], sp, ddf[ia]); + spline(p, f1[ia], sp, ddf1[ia]); + spline(p, f_th[ia], sp, ddf_th[ia]); + spline(p, f1_th[ia], sp, ddf1_th[ia]); + for (ip = 0; ip < sp; ip++) + { + ddf_ud[ia][sp - ip - 1] = ddf[ia][ip]; + ddf1_ud[ia][sp - ip - 1] = ddf1[ia][ip]; + ddf_th_ud[ia][sp - ip - 1] = ddf_th[ia][ip]; + ddf1_th_ud[ia][sp - ip - 1] = ddf1_th[ia][ip]; + } + } + + for (ia = 0; ia < sa; ia++) + { + termC[ia] = (ddf[ia][0] * (3 * p[1] - p[0]) + ddf[ia][sp - 1] * (p[sp - 1] - 3.0 * p[sp - 2])) / 4.0; + for (ip = 0; ip < sp; ip++) + { + termC[ia] += dp6 * ddf[ia][ip]; + } + termC_th[ia] = (ddf_th[ia][0] * (3 * p[1] - p[0]) + ddf_th[ia][sp - 1] * (p[sp - 1] - 3.0 * p[sp - 2])) / 4.0; + for (ip = 0; ip < sp; ip++) + { + termC_th[ia] += dp6 * ddf_th[ia][ip]; + } + } + + // Starting the calculation of ff(x1,x2). + for (ix1 = 0; ix1 < sx2; ix1++) + { + for (ix2 = 0; ix2 <= ix1; ix2++) + { + aux = sqrt(1.0 - x2[ix2] * x2[ix2]); + if (fabs(x2[ix2]) >= 1.0 || fabs(x1[ix1]) >= aux) + { + continue; + } + + // Computation of h_rho + x = x2[ix2] * cos(th[ith]) - x1[ix1] * sin(th[ith]); + for (ip = 0; ip < sp; ip++) + { + double val = fabs(x - p[ip]); + lg[ip] = val < 2e-6 ? 0. : std::log(val); // Using std::log to specify the namespace + } + + for (ia = 0; ia < sa; ia++) + { + image_pos = ia; // 2*ia; + + image[image_pos][image_min_x + sx - ix1 - 1][image_min_y + ix2] + += hilbert_der(x, f[ia], ddf[ia], p, sp, lg, termC[ia]) * c_int; // bot-left + if (ix2 < sy2 - 1) + { + image[image_pos][image_min_x + sx - ix1 - 1][image_min_y + sy - ix2 - 1] + += hilbert_der(x, f1[ia], ddf1[ia], p, sp, lg, termC[ia]) * c_int; // bot-right + } + if (ix1 < sx2 - 1) + { + image[image_pos][image_min_x + ix1][image_min_y + ix2] + -= hilbert_der(-x, f1_ud[ia], ddf1_ud[ia], p_ud, sp, lg, termC[ia]) * c_int; // top-left + } + if ((ix1 < sx2 - 1) && (ix2 < sy2 - 1)) + { + image[image_pos][image_min_x + ix1][image_min_y + sy - ix2 - 1] + -= hilbert_der(-x, f_ud[ia], ddf_ud[ia], p_ud, sp, lg, termC[ia]) * c_int; // top-right + } + + if (ith <= sth2 && ix1 != ix2) + { + image[image_pos][image_min_x + sx - ix2 - 1][image_min_y + ix1] + -= hilbert_der(-x, f_th_ud[ia], ddf_th_ud[ia], p_ud, sp, lg, termC_th[ia]) * c_int; // bot-left + if (ix2 < sy2 - 1) + { + image[image_pos][image_min_x + ix2][image_min_y + ix1] + += hilbert_der(x, f1_th[ia], ddf1_th[ia], p, sp, lg, termC_th[ia]) * c_int; // bot-right + } + if (ix1 < sx2 - 1) + { + image[image_pos][image_min_x + sx - ix2 - 1][image_min_y + sx - ix1 - 1] + -= hilbert_der(-x, f1_th_ud[ia], ddf1_th_ud[ia], p_ud, sp, lg, termC_th[ia]) * c_int; // top-left + } + if ((ix1 < sx2 - 1) && (ix2 < sy2 - 1)) + { + image[image_pos][image_min_x + ix2][image_min_y + sx - ix1 - 1] + += hilbert_der(x, f_th[ia], ddf_th[ia], p, sp, lg, termC_th[ia]) * c_int; // top-right + } + } + else if (ith > sth2 && ix1 != ix2) + { + image[image_pos][image_min_x + sx - ix2 - 1][image_min_y + ix1] + += hilbert_der(x, f_th[ia], ddf_th[ia], p, sp, lg, termC_th[ia]) * c_int; // bot-left + if (ix2 < sy2 - 1) + { + image[image_pos][image_min_x + ix2][image_min_y + ix1] + -= hilbert_der(-x, f1_th_ud[ia], ddf1_th_ud[ia], p_ud, sp, lg, termC_th[ia]) * c_int; // bot-right + } + if (ix1 < sx2 - 1) + { + image[image_pos][image_min_x + sx - ix2 - 1][image_min_y + sx - ix1 - 1] + += hilbert_der(x, f1_th[ia], ddf1_th[ia], p, sp, lg, termC_th[ia]) * c_int; // top-left + } + if ((ix1 < sx2 - 1) && (ix2 < sy2 - 1)) + { + image[image_pos][image_min_x + ix2][image_min_y + sx - ix1 - 1] + -= hilbert_der(-x, f_th_ud[ia], ddf_th_ud[ia], p_ud, sp, lg, termC_th[ia]) * c_int; // top-right + } + } + } + } + } + } + + return Succeeded::yes; +} + +float +SRT2DReconstruction::hilbert_der(float x, + const std::vector& f, + const std::vector& ddf, + const std::vector& p, + int sp, + const std::vector& lg, + float termC) const +{ + float term, trm0, termd0; + + const float d = p[1] - p[0]; + const float d_div_6 = d / 6.0; + const float minus_half_div_d = -0.5 / d; + + term = 0.5 * (ddf[sp - 2] - ddf[0]) * x + termC; + + term += ((f[sp - 1] - f[sp - 2]) / d + ddf[sp - 2] * (d_div_6 + minus_half_div_d * (p[sp - 1] - x) * (p[sp - 1] - x)) + + ddf[sp - 1] * (-d_div_6 - minus_half_div_d * (p[sp - 2] - x) * (p[sp - 2] - x))) + * lg[sp - 1]; + + trm0 = d_div_6 + minus_half_div_d * (p[1] - x) * (p[1] - x); + termd0 = (f[1] - f[0]) / d + ddf[0] * trm0 + ddf[1] * (-d_div_6 - minus_half_div_d * (p[0] - x) * (p[0] - x)); + + term -= termd0 * lg[0]; + + for (int ip = 0; ip < sp - 2; ip++) + { + float trm1 = d_div_6 + minus_half_div_d * (p[ip + 2] - x) * (p[ip + 2] - x); + float termd = (f[ip + 2] - f[ip + 1]) / d + ddf[ip + 1] * trm1 - ddf[ip + 2] * trm0; + term += (termd0 - termd) * lg[ip + 1]; + termd0 = termd; + trm0 = trm1; + } + + return term; +} + +void +SRT2DReconstruction::spline(const std::vector& x, const std::vector& y, int n, std::vector& y2) const +{ + int i, k; + float qn, un; + std::vector u(n); + y2[0] = 0.0; + u[0] = 0.0; + for (i = 1; i < n - 1; i++) + { + float sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]); + float p = sig * y2[i - 1] + 2.0; + y2[i] = (sig - 1.0) / p; + u[i] = (6.0 * ((y[i + 1] - y[i]) / (x[i + 1] - x[i]) - (y[i] - y[i - 1]) / (x[i] - x[i - 1])) / (x[i + 1] - x[i - 1]) + - sig * u[i - 1]) + / p; + } + qn = 0.0; + un = 0.0; + y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0); + for (k = n - 2; k >= 0; k--) + y2[k] = y2[k] * y2[k + 1] + u[k]; + return; +} + +float +SRT2DReconstruction::integ(float dist, int max, const std::vector& ff) const +{ + int k, intg; + intg = ff[0]; + for (k = 1; k < max; k++) + { + intg += ff[k]; + } + return intg * dist / max; +} + +END_NAMESPACE_STIR diff --git a/src/analytic/SRT2DSPECT/CMakeLists.txt b/src/analytic/SRT2DSPECT/CMakeLists.txt new file mode 100644 index 0000000000..e09288b330 --- /dev/null +++ b/src/analytic/SRT2DSPECT/CMakeLists.txt @@ -0,0 +1,23 @@ +# +# + +set(dir analytic_SRT2DSPECT) + +set (dir_LIB_SOURCES ${dir}_LIB_SOURCES) + +set(${dir_LIB_SOURCES} + SRT2DSPECTReconstruction.cxx +) + +#$(dir)_REGISTRY_SOURCES:= + +include(stir_lib_target) +target_link_libraries(analytic_SRT2DSPECT recon_buildblock IO ) + +set (dir_EXE_SOURCES ${dir}_EXE_SOURCES) + +set(${dir_EXE_SOURCES} + SRT2DSPECT.cxx +) + +include(stir_exe_targets) diff --git a/src/analytic/SRT2DSPECT/SRT2DSPECT.cxx b/src/analytic/SRT2DSPECT/SRT2DSPECT.cxx new file mode 100644 index 0000000000..c73e1cc3bd --- /dev/null +++ b/src/analytic/SRT2DSPECT/SRT2DSPECT.cxx @@ -0,0 +1,36 @@ +// +// +/*! + \file + \ingroup analytic + \brief Main program for SRT2DSPECT reconstruction + \author Dimitra Kyriakopoulou + \author Kris Thielemans +*/ +/* + Copyright (C) 2024, University College London + + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +#include "stir/analytic/SRT2DSPECT/SRT2DSPECTReconstruction.h" +#include "stir/Succeeded.h" +#ifndef PARALLEL +# define Main main +#else +# define Main master_main +#endif + +USING_NAMESPACE_STIR + +int +Main(int argc, char** argv) +{ + SRT2DSPECTReconstruction reconstruction_object(argc > 1 ? argv[1] : ""); + + return reconstruction_object.reconstruct() == Succeeded::yes ? EXIT_SUCCESS : EXIT_FAILURE; +} diff --git a/src/analytic/SRT2DSPECT/SRT2DSPECTReconstruction.cxx b/src/analytic/SRT2DSPECT/SRT2DSPECTReconstruction.cxx new file mode 100644 index 0000000000..ecf3874dca --- /dev/null +++ b/src/analytic/SRT2DSPECT/SRT2DSPECTReconstruction.cxx @@ -0,0 +1,761 @@ +/* + Copyright (C) 2024 University College London + + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +/*! + \file + \ingroup analytic + \brief Implementation of class stir::SRT2DSPECTReconstruction + + \author Dimitra Kyriakopoulou + \author Kris Thielemans +*/ + +#include "stir/analytic/SRT2DSPECT/SRT2DSPECTReconstruction.h" +#include "stir/VoxelsOnCartesianGrid.h" +#include "stir/ProjDataInfoCylindricalArcCorr.h" +#include "stir/SSRB.h" +#include "stir/ProjDataInMemory.h" +#include "stir/Array.h" +#include +#include "stir/Sinogram.h" +#include "stir/Viewgram.h" +#include "stir/Bin.h" +#include "stir/round.h" +#include "stir/display.h" +#include +#include "stir/IO/interfile.h" +#include "stir/info.h" +#include + +#include "stir/SegmentByView.h" +#include "stir/ArcCorrection.h" +#include "stir/shared_ptr.h" + +/*#ifdef STIR_OPENMP +# include +#endif*/ + +#include // For M_PI and other math functions +#ifndef M_PI +# define M_PI 3.14159265358979323846 +#endif + +#ifdef STIR_OPENMP +# include "stir/num_threads.h" +#endif + +#include "stir/Coordinate3D.h" + +START_NAMESPACE_STIR + +const char* const SRT2DSPECTReconstruction::registered_name = "SRT2DSPECT"; + +void +SRT2DSPECTReconstruction::set_defaults() +{ + base_type::set_defaults(); + attenuation_projection_filename = ""; + num_segments_to_combine = -1; +} + +void +SRT2DSPECTReconstruction::initialise_keymap() +{ + base_type::initialise_keymap(); + + parser.add_start_key("SRT2DSPECTParameters"); + parser.add_stop_key("End"); + parser.add_key("num_segments_to_combine with SSRB", &num_segments_to_combine); + parser.add_key("attenuation projection filename", &attenuation_projection_filename); +} + +void +SRT2DSPECTReconstruction::ask_parameters() +{ + base_type::ask_parameters(); + num_segments_to_combine = ask_num("num_segments_to_combine (must be odd)", -1, 101, -1); + attenuation_projection_filename = ask_string("attenuation projection filename"); +} + +bool +SRT2DSPECTReconstruction::post_processing() +{ + return base_type::post_processing(); +} + +Succeeded +SRT2DSPECTReconstruction::set_up(shared_ptr const& target_data_sptr) +{ + if (base_type::set_up(target_data_sptr) == Succeeded::no) + return Succeeded::no; + atten_data_ptr = ProjData::read_from_file(attenuation_projection_filename); + + if (num_segments_to_combine >= 0 && num_segments_to_combine % 2 == 0) + error(boost::format("num_segments_to_combine has to be odd (or -1), but is %d") % num_segments_to_combine); + + if (num_segments_to_combine == -1) + { + const shared_ptr proj_data_info_cyl_sptr + = dynamic_pointer_cast(proj_data_ptr->get_proj_data_info_sptr()); + + if (is_null_ptr(proj_data_info_cyl_sptr)) + num_segments_to_combine = 1; // cannot SSRB non-cylindrical data yet + else + { + if (proj_data_info_cyl_sptr->get_min_ring_difference(0) != proj_data_info_cyl_sptr->get_max_ring_difference(0) + || proj_data_info_cyl_sptr->get_num_segments() == 1) + num_segments_to_combine = 1; + else + num_segments_to_combine = 3; + } + } + + return Succeeded::yes; +} + +std::string +SRT2DSPECTReconstruction::method_info() const +{ + return "SRT2DSPECT"; +} + +SRT2DSPECTReconstruction::SRT2DSPECTReconstruction(const std::string& parameter_filename) +{ + initialise(parameter_filename); + info(boost::format("%1%") % parameter_info()); +} + +SRT2DSPECTReconstruction::SRT2DSPECTReconstruction() +{ + set_defaults(); +} + +SRT2DSPECTReconstruction::SRT2DSPECTReconstruction(const shared_ptr& proj_data_ptr_v, + const int num_segments_to_combine_v) +{ + set_defaults(); + proj_data_ptr = proj_data_ptr_v; + num_segments_to_combine = num_segments_to_combine_v; +} + +Succeeded +SRT2DSPECTReconstruction::actual_reconstruct(shared_ptr> const& density_ptr) +{ + + // perform SSRB + if (num_segments_to_combine > 1) + { + const ProjDataInfoCylindrical& proj_data_info_cyl + = dynamic_cast(*proj_data_ptr->get_proj_data_info_sptr()); + + // full_log << "SSRB combining " << num_segments_to_combine + // << " segments in input file to a new segment 0\n" << std::endl; + + shared_ptr ssrb_info_sptr( + SSRB(proj_data_info_cyl, num_segments_to_combine, 1, 0, (num_segments_to_combine - 1) / 2)); + shared_ptr proj_data_to_SRT_ptr(new ProjDataInMemory(proj_data_ptr->get_exam_info_sptr(), ssrb_info_sptr)); + SSRB(*proj_data_to_SRT_ptr, *proj_data_ptr); + proj_data_ptr = proj_data_to_SRT_ptr; + } + else + { + // just use the proj_data_ptr we have already + } + + // check if segment 0 has direct sinograms + { + const float tan_theta = proj_data_ptr->get_proj_data_info_sptr()->get_tantheta(Bin(0, 0, 0, 0)); + if (fabs(tan_theta) > 1.E-4) + { + warning("SRT2DSPECT: segment 0 has non-zero tan(theta) %g", tan_theta); + return Succeeded::no; + } + } + + auto pdi_sptr = dynamic_pointer_cast(proj_data_ptr->get_proj_data_info_sptr()); + if (!pdi_sptr) + { + error("SPECT data should correspond to ProjDataInfoCylindricalArcCorr"); + } + + VoxelsOnCartesianGrid& image = dynamic_cast&>(*density_ptr); + density_ptr->fill(0); + Sinogram sino = proj_data_ptr->get_empty_sinogram(0, 0); + // Viewgram view = proj_data_ptr->get_empty_viewgram(0, 0); + Viewgram view = proj_data_ptr->get_viewgram(0, 0); + + Viewgram view_atten = atten_data_ptr->get_empty_viewgram(0, 0); + + // Retrieve runtime-dependent sizes + const int sp = view.get_num_tangential_poss(); // const int sp = proj_data_ptr->get_num_tangential_poss(); + const int sth = proj_data_ptr->get_num_views(); + const int sa = proj_data_ptr->get_num_axial_poss(0); + + // RelatedViewgrams viewgrams; + + const int sx = image.get_x_size(); + const int sy = image.get_y_size(); + + // c ---------------------------------------------- + // c The rest of the variables used by the program. + // c ---------------------------------------------- + int i, j, k1, k2; + int ith, ia, ip, ix1, ix2; // extra + float aux, a, b, f_node; + float x; // extra + + const int image_min_x = image.get_min_x(); + const int image_min_y = image.get_min_y(); + + // Projection angles and sampling arrays + std::vector th(sth, 0); // Theta values for each view angle. + std::vector p(sp, 0); // Tangential positions for projections. + std::vector x1(sx, 0), x2(sy, 0); // Normalized coordinates for image reconstruction. + + // Projection data storage + std::vector> g(sa, std::vector(sp, 0)); // Projection data matrix. + std::vector> ddg(sa, + std::vector(sp, 0)); // Second derivative of projections for spline interpolation. + + // Number of angular divisions + const int Nt = 8; // Number of angular subdivisions used for integration. + const int Nmul = sth / Nt; // Multiplier for view angle processing. + + std::vector lg(sp, 0); // Logarithmic values for tangential positions. + float dh1[Nt], dh2[Nt]; // Derivative values at angular subdivisions. + float t[Nt]; // Angular subdivision points. + + // Intermediate Hilbert transform results for each axial position + std::vector> hilb(sa, std::vector(sp, 0)); // Hilbert transform results. + std::vector> fcpe(sa, std::vector(sp, 0)); // Cosine of phase exponentials. + std::vector> fspe(sa, std::vector(sp, 0)); // Sine of phase exponentials. + std::vector> fc(sa, std::vector(sp, 0)); // Cosine-filtered projections. + std::vector> fs(sa, std::vector(sp, 0)); // Sine-filtered projections. + std::vector> ddfc(sa, std::vector(sp, 0)); // Second derivative for cosine-filtered projections. + std::vector> ddfs(sa, std::vector(sp, 0)); // Second derivative for sine-filtered projections. + + // Storage for second derivatives and interpolations + std::vector> f(sa, std::vector(sp, 0)); // Attenuation projections. + std::vector> ddf(sa, std::vector(sp, 0)); // Second derivatives of attenuation projections. + + // Variables for Hilbert transform and interpolation results + float rho, h, fcme_fin, fsme_fin, fc_fin, fs_fin, fcpe_fin, fspe_fin, hc_fin, + hs_fin; // Variables for Hilbert transform calculations. + float I, Ft1, Ft2, rho1, rho2, tau, tau1, tau2, rx1, rx2; // Variables for intermediate calculations. + float gx, w, F; // Variables for integration and filtering calculations. + + // Cache for logarithmic differences used in Hilbert transforms + std::vector> lg1_cache(Nt / 2, std::vector(sp - 1, 0)); // Logarithmic differences for \a rho1. + std::vector> lg2_cache(Nt / 2, std::vector(sp - 1, 0)); // Logarithmic differences for \a rho2. + + // 3D array for storing reconstructed image slices + IndexRange<3> range(Coordinate3D(0, 0, 0), Coordinate3D(sa - 1, sx - 1, sy - 1)); + Array<3, float> rx1x2th(range); + rx1x2th.fill(0.0F); // Initialize the reconstruction array to zeros. + + // Caches for Hilbert transforms in multiple angles + std::vector>> f_cache( + sa, std::vector>(Nt / 2, std::vector(sp, 0))); // Cache for filtered projections. + std::vector>> ddf_cache( + sa, std::vector>(Nt / 2, std::vector(sp, 0))); // Cache for second derivatives. + std::vector>> f1_cache( + sa, std::vector>(Nt / 2, std::vector(sp, 0))); // Cache for mirrored projections. + std::vector>> ddf1_cache( + sa, + std::vector>(Nt / 2, + std::vector(sp, 0))); // Cache for second derivatives of mirrored projections. + + /* #ifdef STIR_OPENMP + set_num_threads(); + #pragma omp single + info("Using OpenMP-version of SRT2D with " + std::to_string(omp_get_num_threads()) + + " threads on " + std::to_string(omp_get_num_procs()) + " processors."); + #endif*/ + + // c -------------------------- + // c Put theta and p in arrays. + // c -------------------------- + for (i = 0; i < sth; i++) + th[i] = i * 2 * M_PI / sth; + for (int it = 0; it < Nt; it++) + t[it] = it * 2 * M_PI / Nt; + for (j = 0; j < sp; j++) + p[j] = -1.0 + 2.0 * j / (sp - 1); + + // c ------------------------ + // c Put x1 and x2 in arrays. + // c ------------------------ + for (k1 = 0; k1 < sx; k1++) + x1[k1] = -1.0 + 2.0 * k1 / (sx - 1); + for (k2 = 0; k2 < sx; k2++) + x2[k2] = -1.0 + 2.0 * k2 / (sx - 1); + + for (int it = 0; it < Nt / 2; it++) + { + view_atten = atten_data_ptr->get_viewgram(Nmul * it, 0); + for (int ia = 0; ia < sa; ia++) + { + for (int ip = 0; ip < sp; ip++) + { + f_cache[ia][it][ip] + = view_atten[view_atten.get_min_axial_pos_num() + ia][view_atten.get_min_tangential_pos_num() + ip]; //*.15; + } + } + for (int ia = 0; ia < sa; ia++) + { + spline(p, f_cache[ia][it], sp, ddf_cache[ia][it]); + } + + for (int ia = 0; ia < sa; ia++) + { + for (int ip = 0; ip < sp; ip++) + { + f1_cache[ia][it][sp - ip - 1] = f_cache[ia][it][ip]; + } + } + for (int ia = 0; ia < sa; ia++) + { + for (int ip = 0; ip < sp; ip++) + { + ddf1_cache[ia][it][sp - ip - 1] = ddf_cache[ia][it][ip]; + } + } + } + + //-- Starting calculations per view + // 2D algorithm only + // At the moment, the parallelization produces artifacts that the non-parallelized version does not have. That's why it's + // commented out. + /*#ifdef STIR_OPENMP + # pragma omp parallel firstprivate(f, ddf, f_cache, ddf_cache, f1_cache, ddf1_cache, hilb, fcpe, fspe, fc, fs, ddfc, ddfs, aux, + rho, lg, tau, a, b, tau1, tau2, w, rho1, rho2, lg1_cache, lg2_cache, f_node, h, fcme_fin, fsme_fin, fcpe_fin, fspe_fin, gx, + fc_fin, fs_fin, hc_fin, hs_fin, dh1, dh2, Ft1, Ft2, F, I, rx1, rx2) \ shared(view, view_atten, do_arc_correction, + arc_correction, p, th, x1, x2, image, proj_data_ptr, atten_data_ptr, rx1x2th) private(ith, ia, ip, ix1, ix2) # pragma omp for + schedule(dynamic) nowait #endif */ + for (ith = 0; ith < sth; ith++) + { + info(boost::format("View %d of %d") % ith % sth); + + //-- Loading the viewgram + /*#ifdef STIR_OPENMP + # pragma omp critical + #endif*/ + { + view = proj_data_ptr->get_viewgram(ith, 0); + view_atten = atten_data_ptr->get_viewgram(ith, 0); + for (ia = 0; ia < sa; ia++) + { + + for (ip = 0; ip < sp; ip++) + { + g[ia][ip] = view[view.get_min_axial_pos_num() + ia][view.get_min_tangential_pos_num() + ip]; + f[ia][ip] + = view_atten[view_atten.get_min_axial_pos_num() + ia][view_atten.get_min_tangential_pos_num() + ip] * 0.1; + } + } + } + //-- Calculation of second derivative by use of function spline + for (ia = 0; ia < sa; ia++) + { + spline(p, g[ia], sp, ddg[ia]); + spline(p, f[ia], sp, ddf[ia]); + } + + //---- calculate h(rho,theta) for all rho, theta + for (ia = 0; ia < sa; ia++) + { + for (ip = 0; ip < sp; ip++) + { + hilb[ia][ip] = hilbert_node(p[ip], f[ia], ddf[ia], p, sp, f[ia][ip]); + + fcpe[ia][ip] = exp(0.5 * f[ia][ip]) * cos(hilb[ia][ip] / (2 * M_PI)); + fspe[ia][ip] = exp(0.5 * f[ia][ip]) * sin(hilb[ia][ip] / (2 * M_PI)); + + fc[ia][ip] = fcpe[ia][ip] * g[ia][ip]; + fs[ia][ip] = fspe[ia][ip] * g[ia][ip]; + } + //-- calculate ddfc, ddfs for all \rho, \theta + spline(p, fc[ia], sp, ddfc[ia]); + spline(p, fs[ia], sp, ddfs[ia]); + } + + //---- calculate r(x1, x2, theta) + for (ix1 = 0; ix1 < sx; ix1++) + { + for (ix2 = 0; ix2 < sy; ix2++) + { + aux = sqrt(1. - x2[ix2] * x2[ix2]); + if (fabs(x2[ix2]) >= 1. || fabs(x1[ix1]) >= aux) + continue; + + rho = x2[ix2] * cos(th[ith]) - x1[ix1] * sin(th[ith]); + + int i = floor((rho + 1) * (sp - 1) / 2); + float p1 = p[i]; + float p2 = p[i + 1]; + float A = (p2 - rho) / (p2 - p1); + float B = 1 - A; + float C = 1.0 / 6 * (A * A * A - A) * (p2 - p1) * (p2 - p1); + float D = 1.0 / 6 * (B * B * B - B) * (p2 - p1) * (p2 - p1); + + for (ip = 0; ip < sp; ip++) + { + double val = fabs(rho - p[ip]); + lg[ip] = val < 2e-6 ? 0. : std::log(val); // Using std::log to specify the namespace + } + + // calculate I + tau = x2[ix2] * sin(th[ith]) + x1[ix1] * cos(th[ith]); + if (tau >= 0) + { + a = tau; + b = sqrt(1 - rho * rho); + } + else + { + a = -sqrt(1 - rho * rho); + b = tau; + } + + tau1 = a + (b - a) * (1.0 / 2 - sqrt(3.0) / 6); + tau2 = a + (b - a) * (1.0 / 2 + sqrt(3.0) / 6); + w = 1.0 / 2 * (b - a); + + for (int it = 0; it < Nt / 2; it++) + { + rho1 = tau1 * sin(th[ith] - t[it]) + rho * cos(th[ith] - t[it]); + rho2 = tau2 * sin(th[ith] - t[it]) + rho * cos(th[ith] - t[it]); + + for (ip = 0; ip < sp - 1; ip++) + { + lg1_cache[it][ip] = log(fabs((p[ip + 1] - rho1) / (p[ip] - rho1))); + if (fabs(p[ip + 1] - rho1) < 2e-6 || fabs(p[ip] - rho1) < 2e-6) + lg1_cache[it][ip] = 0.; + lg2_cache[it][ip] = log(fabs((p[ip + 1] - rho2) / (p[ip] - rho2))); + if (fabs(p[ip + 1] - rho2) < 2e-6 || fabs(p[ip] - rho2) < 2e-6) + lg2_cache[it][ip] = 0.; + } + } + + for (ia = 0; ia < sa; ia++) + { + // Example of how to choose particular slices to be reconstructed + // if(ia!=20 && ia!=31 && ia!=70 && ia!=71 &&ia!=81 && ia!=100) continue; + + f_node = A * f[ia][i] + B * f[ia][i + 1] + C * ddf[ia][i] + D * ddf[ia][i + 1]; + + // calculate fcme, fsme, fc, fs, hc, hs + + h = hilbert(rho, f[ia], ddf[ia], p, sp, lg); + fcme_fin = exp(-0.5 * f_node) * cos(h / (2 * M_PI)); + fsme_fin = exp(-0.5 * f_node) * sin(h / (2 * M_PI)); + + fcpe_fin = exp(0.5 * f_node) * cos(h / (2 * M_PI)); + fspe_fin = exp(0.5 * f_node) * sin(h / (2 * M_PI)); + + gx = splint(p, g[ia], ddg[ia], sp, rho); + + fc_fin = fcpe_fin * gx; + fs_fin = fspe_fin * gx; + + hc_fin = hilbert(rho, fc[ia], ddfc[ia], p, sp, lg); + hs_fin = hilbert(rho, fs[ia], ddfs[ia], p, sp, lg); + + rx1x2th[ia][ix1][ix2] = fcme_fin * (1.0 / M_PI * hc_fin + fs_fin) + fsme_fin * (1.0 / M_PI * hs_fin - fc_fin); + + // calculate I + for (int it = 0; it < Nt / 2; it++) + { + rho1 = tau1 * sin(th[ith] - t[it]) + rho * cos(th[ith] - t[it]); + rho2 = tau2 * sin(th[ith] - t[it]) + rho * cos(th[ith] - t[it]); + hilbert_der_double(rho1, + f_cache[ia][it], + ddf_cache[ia][it], + f1_cache[ia][it], + ddf1_cache[ia][it], + p, + sp, + &dh1[it], + &dh1[it + Nt / 2], + lg1_cache[it]); + hilbert_der_double(rho2, + f_cache[ia][it], + ddf_cache[ia][it], + f1_cache[ia][it], + ddf1_cache[ia][it], + p, + sp, + &dh2[it], + &dh2[it + Nt / 2], + lg2_cache[it]); + } + + Ft1 = -1.0 / (4.0 * M_PI * M_PI) * integ(2 * M_PI, Nt, dh1); + Ft2 = -1.0 / (4.0 * M_PI * M_PI) * integ(2 * M_PI, Nt, dh2); + F = w * Ft1 + w * Ft2; + + I = exp(f_node - F); + + rx1x2th[ia][ix1][ix2] = I * rx1x2th[ia][ix1][ix2]; + } + } + } + + //---- calculate g(x1, x2) + for (ia = 0; ia < sa; ia++) + { + for (ix1 = 0; ix1 < sx; ix1++) + { + for (ix2 = 0; ix2 < sy; ix2++) + { + aux = sqrt(1.0 - x2[ix2] * x2[ix2]); + if (fabs(x2[ix2]) >= 1.0 || fabs(x1[ix1]) >= aux) + { + continue; + } + + if (x1[ix1] < 0) + { + rx1 = (-3.0 * rx1x2th[ia][ix1][ix2] + 4.0 * rx1x2th[ia][ix1 + 1][ix2] - rx1x2th[ia][ix1 + 2][ix2]) + / (2.0 * (2.0 / (sx - 1))); + } + else + { + rx1 = (3.0 * rx1x2th[ia][ix1][ix2] - 4.0 * rx1x2th[ia][ix1 - 1][ix2] + rx1x2th[ia][ix1 - 2][ix2]) + / (2.0 * (2.0 / (sx - 1))); + } + + if (x2[ix2] < 0) + { + rx2 = (-3.0 * rx1x2th[ia][ix1][ix2] + 4.0 * rx1x2th[ia][ix1][ix2 + 1] - rx1x2th[ia][ix1][ix2 + 2]) + / (2.0 * (2.0 / (sy - 1))); + } + else + { + rx2 = (3.0 * rx1x2th[ia][ix1][ix2] - 4.0 * rx1x2th[ia][ix1][ix2 - 1] + rx1x2th[ia][ix1][ix2 - 2]) + / (2.0 * (2.0 / (sy - 1))); + } + + /*#ifdef STIR_OPENMP + # pragma omp critical + #endif*/ + { + image[ia][image_min_x + sx - ix1 - 1][image_min_y + ix2] + += 1.0 / (4.0 * M_PI) * (rx1 * sin(th[ith]) - rx2 * cos(th[ith])) * (2.0 * M_PI / sth) * 6.23; + } + } + } + } + } // slice + + return Succeeded::yes; +} + +float +SRT2DSPECTReconstruction::hilbert_node( + float x, const std::vector& f, const std::vector& ddf, const std::vector& p, int sp, float fn) const +{ + float dh; + + dh = 0; + for (int i = 0; i < sp - 1; i++) + { + dh = dh - f[i] + f[i + 1] + + 1.0 / 36 + * (4 * p[i] * p[i] - 5 * p[i] * p[i + 1] - 5 * p[i + 1] * p[i + 1] - 3 * (p[i] - 5 * p[i + 1]) * x - 6 * x * x) + * ddf[i] + + 1.0 / 36 + * (5 * p[i] * p[i] + 5 * p[i] * p[i + 1] - 4 * p[i + 1] * p[i + 1] - 3 * (5 * p[i] - p[i + 1]) * x + 6 * x * x) + * ddf[i + 1]; + } + + if (fabs(x) == 1) + { + dh = 2 / (sp - 1) * dh; + } + else + { + dh = fn * log((1 - x) / (1 + x)) + 2 / (sp - 1) * dh; + } + + return dh; +} + +float +SRT2DSPECTReconstruction::hilbert(float x, + const std::vector& f, + const std::vector& ddf, + const std::vector& p, + int sp, + std::vector& lg) const +{ + float dh, Di; + int i; + + i = 0; + Di = -1.0 / (p[i] - p[i + 1]) + * ((p[i + 1] - x) * f[i] - (p[i] - x) * f[i + 1] + - 1.0 / 6 * (p[i] - x) * (p[i + 1] - x) + * ((p[i] - 2 * p[i + 1] + x) * ddf[i] + (2 * p[i] - p[i + 1] - x) * ddf[i + 1])); + dh = -f[i] + f[i + 1] + + 1.0 / 36 * (4 * p[i] * p[i] - 5 * p[i] * p[i + 1] - 5 * p[i + 1] * p[i + 1] - 3 * (p[i] - 5 * p[i + 1]) * x - 6 * x * x) + * ddf[i] + + 1.0 / 36 * (5 * p[i] * p[i] + 5 * p[i] * p[i + 1] - 4 * p[i + 1] * p[i + 1] - 3 * (5 * p[i] - p[i + 1]) * x + 6 * x * x) + * ddf[i + 1] + - Di * lg[i]; + + for (i = 1; i < sp - 2; i++) + { + + float Di1 = -1.0 / (p[i] - p[i + 1]) + * ((p[i + 1] - x) * f[i] - (p[i] - x) * f[i + 1] + - 1.0 / 6 * (p[i] - x) * (p[i + 1] - x) + * ((p[i] - 2 * p[i + 1] + x) * ddf[i] + (2 * p[i] - p[i + 1] - x) * ddf[i + 1])); + + dh = dh - f[i] + f[i + 1] + + 1.0 / 36 + * (4 * p[i] * p[i] - 5 * p[i] * p[i + 1] - 5 * p[i + 1] * p[i + 1] - 3 * (p[i] - 5 * p[i + 1]) * x - 6 * x * x) + * ddf[i] + + 1.0 / 36 + * (5 * p[i] * p[i] + 5 * p[i] * p[i + 1] - 4 * p[i + 1] * p[i + 1] - 3 * (5 * p[i] - p[i + 1]) * x + 6 * x * x) + * ddf[i + 1] + + (Di - Di1) * lg[i + 1]; + + Di = Di1; + } + + i = sp - 2; + Di = -1.0 / (p[i] - p[i + 1]) + * ((p[i + 1] - x) * f[i] - (p[i] - x) * f[i + 1] + - 1.0 / 6 * (p[i] - x) * (p[i + 1] - x) + * ((p[i] - 2 * p[i + 1] + x) * ddf[i] + (2 * p[i] - p[i + 1] - x) * ddf[i + 1])); + dh = dh - f[i] + f[i + 1] + + 1.0 / 36 * (4 * p[i] * p[i] - 5 * p[i] * p[i + 1] - 5 * p[i + 1] * p[i + 1] - 3 * (p[i] - 5 * p[i + 1]) * x - 6 * x * x) + * ddf[i] + + 1.0 / 36 * (5 * p[i] * p[i] + 5 * p[i] * p[i + 1] - 4 * p[i + 1] * p[i + 1] - 3 * (5 * p[i] - p[i + 1]) * x + 6 * x * x) + * ddf[i + 1] + + Di * lg[sp - 1]; + + dh = 2.0 / (sp - 1) * dh; + + return dh; +} + +void +SRT2DSPECTReconstruction::hilbert_der_double(float x, + const std::vector& f, + const std::vector& ddf, + const std::vector& f1, + const std::vector& ddf1, + const std::vector& p, + int sp, + float* dhp, + float* dh1p, + const std::vector& lg) const +{ + float dh, dh1, dp; + dh = 0; + dh1 = 0; + dp = p[1] - p[2]; // float pix = 0., pi1x = 0.; + for (int i = 0; i < sp - 1; i++) + { + float pix = fabs(p[i] - x) > 2e-6 ? f[i] / (p[i] - x) : 0.; + float pi1x = fabs(p[i + 1] - x) > 2e-6 ? f[i + 1] / (p[i + 1] - x) : 0.; + dh = dh + pix - pi1x - 1.0 / 4 * (p[i] - 3 * p[i + 1] + 2 * x) * ddf[i] + - 1.0 / 4 * (3 * p[i] - p[i + 1] - 2 * x) * ddf[i + 1] + + ((f[i] - f[i + 1]) / dp - 1.0 / 6 * (p[i] - p[i + 1] - (3 * (p[i + 1] - x) * (p[i + 1] - x)) / dp) * ddf[i] + + 1.0 / 6 * (p[i] - p[i + 1] - (3 * (p[i] - x) * (p[i] - x)) / dp) * ddf[i + 1]) + * lg[i]; + pix = fabs(p[i] - x) > 2e-6 ? f1[i] / (p[i] - x) : 0.; + pi1x = fabs(p[i + 1] - x) > 2e-6 ? f1[i + 1] / (p[i + 1] - x) : 0.; + dh1 = dh1 + pix - pi1x - 1.0 / 4 * (p[i] - 3 * p[i + 1] + 2 * x) * ddf1[i] + - 1.0 / 4 * (3 * p[i] - p[i + 1] - 2 * x) * ddf1[i + 1] + + ((f1[i] - f1[i + 1]) / dp - 1.0 / 6 * (p[i] - p[i + 1] - (3 * (p[i + 1] - x) * (p[i + 1] - x)) / dp) * ddf1[i] + + 1.0 / 6 * (p[i] - p[i + 1] - (3 * (p[i] - x) * (p[i] - x)) / dp) * ddf1[i + 1]) + * lg[i]; + } + dh = 2.0 / (sp - 1) * dh; + dh1 = 2.0 / (sp - 1) * dh1; + *dhp = dh; + *dh1p = dh1; +} + +float +SRT2DSPECTReconstruction::splint( + const std::vector& xa, const std::vector& ya, const std::vector& y2a, int n, float x) const +{ + int klo, khi; + float h, a, b, y; + + klo = 1; + khi = n; + while (khi - klo > 1) + { + int k = floor((khi + klo) / 2.0); + if (xa[k] > x) + { + khi = k; + } + else + { + klo = k; + } + } + + h = xa[khi] - xa[klo]; + /* if(h == 0) { + error('bad xa input in splint'); + } */ + a = (xa[khi] - x) / h; + b = (x - xa[klo]) / h; + y = a * ya[klo] + b * ya[khi] + ((a * a * a - a) * y2a[klo] + (b * b * b - b) * y2a[khi]) * (h * h) / 6.0; + + return y; +} + +void +SRT2DSPECTReconstruction::spline(const std::vector& x, const std::vector& y, int n, std::vector& y2) const +{ + // function for nanural qubic spline. + int i, k; + float qn, un; + std::vector u(n); + y2[0] = 0.0; + u[0] = 0.0; + for (i = 1; i < n - 1; i++) + { + float sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]); + float p = sig * y2[i - 1] + 2.0; + y2[i] = (sig - 1.0) / p; + u[i] = (6.0 * ((y[i + 1] - y[i]) / (x[i + 1] - x[i]) - (y[i] - y[i - 1]) / (x[i] - x[i - 1])) / (x[i + 1] - x[i - 1]) + - sig * u[i - 1]) + / p; + } + qn = 0.0; + un = 0.0; + y2[n - 1] = (un - qn * u[n - 2]) / (qn * y2[n - 2] + 1.0); + for (k = n - 2; k >= 0; k--) + y2[k] = y2[k] * y2[k + 1] + u[k]; + return; +} + +float +SRT2DSPECTReconstruction::integ(float dist, int max, float ff[]) const +{ + int k, intg; + intg = ff[0]; + for (k = 1; k < max; k++) + { + intg += ff[k]; + } + return intg * dist / max; +} + +END_NAMESPACE_STIR diff --git a/src/buildblock/CMakeLists.txt b/src/buildblock/CMakeLists.txt index add1c45ebc..03e0f85291 100644 --- a/src/buildblock/CMakeLists.txt +++ b/src/buildblock/CMakeLists.txt @@ -59,6 +59,10 @@ set(${dir_LIB_SOURCES} SeparableGaussianArrayFilter.cxx MedianArrayFilter3D.cxx MedianImageFilter3D.cxx + WienerArrayFilter2D.cxx + WienerImageFilter2D.cxx + GammaArrayFilter2D.cxx + GammaImageFilter2D.cxx MinimalArrayFilter3D.cxx MinimalImageFilter3D.cxx SeparableCartesianMetzImageFilter.cxx diff --git a/src/buildblock/GammaArrayFilter2D.cxx b/src/buildblock/GammaArrayFilter2D.cxx new file mode 100644 index 0000000000..8e8d785795 --- /dev/null +++ b/src/buildblock/GammaArrayFilter2D.cxx @@ -0,0 +1,103 @@ +// +// +/*! + \file + \ingroup Array + \brief Implementations for class stir::GammaArrayFilter2D + + \author Dimitra Kyriakopoulou + \author Kris Thielemans +*/ +/* + Copyright (C) 2024, University College London + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +#include "stir/GammaArrayFilter2D.h" +#include +#include + +START_NAMESPACE_STIR + +template +GammaArrayFilter2D::GammaArrayFilter2D() +{} + +template +void +GammaArrayFilter2D::do_it(Array<3, elemT>& out_array, const Array<3, elemT>& in_array) const +{ + assert(out_array.get_index_range() == in_array.get_index_range()); + + const int sx = out_array[0][0].get_length(); + const int sy = out_array[0].get_length(); + const int sa = out_array.get_length(); + const int min_x = in_array[0][0].get_min_index(); + const int min_y = in_array[0].get_min_index(); + const float targetAverage = 0.25; + + for (int ia = 0; ia < sa; ia++) + { + // Find min and max values in the current slice + float min_val = INFINITY, max_val = -INFINITY; + for (int i = 0; i < sx; i++) + { + for (int j = 0; j < sy; j++) + { + min_val = std::min(in_array[ia][min_x + i][min_y + j], min_val); + max_val = std::max(in_array[ia][min_x + i][min_y + j], max_val); + } + } + + // Normalize the slice + for (int i = 0; i < sx; i++) + for (int j = 0; j < sy; j++) + out_array[ia][min_x + i][min_y + j] = (in_array[ia][min_x + i][min_y + j] - min_val) / (max_val - min_val); + + // Calculate average pixel value + float averagePixelValue = 0.0; + int count = 0; + for (int i = 0; i < sx; i++) + { + for (int j = 0; j < sy; j++) + { + if (std::abs(out_array[ia][min_x + i][min_y + j]) > 0.1) + { + count++; + averagePixelValue += out_array[ia][min_x + i][min_y + j]; + } + } + } + averagePixelValue /= count; + + // Apply gamma correction + float gamma_val = 1.0; + if (averagePixelValue > 0.0) + gamma_val = std::log(targetAverage) / std::log(averagePixelValue); + + for (int i = 0; i < sx; i++) + for (int j = 0; j < sy; j++) + out_array[ia][min_x + i][min_y + j] = std::pow(out_array[ia][min_x + i][min_y + j], gamma_val); + + // Restore original scale + for (int i = 0; i < sx; i++) + for (int j = 0; j < sy; j++) + out_array[ia][min_x + i][min_y + j] = out_array[ia][min_x + i][min_y + j] * (max_val - min_val) + min_val; + } +} + +template +bool +GammaArrayFilter2D::is_trivial() const +{ + return false; +} + +// Explicit template instantiation +template class GammaArrayFilter2D; + +END_NAMESPACE_STIR diff --git a/src/buildblock/GammaImageFilter2D.cxx b/src/buildblock/GammaImageFilter2D.cxx new file mode 100644 index 0000000000..de5704fd76 --- /dev/null +++ b/src/buildblock/GammaImageFilter2D.cxx @@ -0,0 +1,82 @@ +// +// +/*! + \file + \ingroup ImageProcessor + \brief Implementations for class stir::GammaImageFilter2D + + \author Dimitra Kyriakopoulou + \author Kris Thielemans +*/ +/* + Copyright (C) 2024, University College London + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +#include "stir/GammaImageFilter2D.h" +#include "stir/DiscretisedDensity.h" +#include + +START_NAMESPACE_STIR + +template +GammaImageFilter2D::GammaImageFilter2D() +{ + std::cout << "Gamma filter initialization" << std::endl; + set_defaults(); +} + +template +Succeeded +GammaImageFilter2D::virtual_set_up(const DiscretisedDensity<3, elemT>& density) +{ + gamma_filter = GammaArrayFilter2D(); + return Succeeded::yes; +} + +template +void +GammaImageFilter2D::virtual_apply(DiscretisedDensity<3, elemT>& density) const +{ + gamma_filter(density); +} + +template +void +GammaImageFilter2D::virtual_apply(DiscretisedDensity<3, elemT>& out_density, + const DiscretisedDensity<3, elemT>& in_density) const +{ + gamma_filter(out_density, in_density); +} + +template +void +GammaImageFilter2D::set_defaults() +{ + base_type::set_defaults(); +} + +template +void +GammaImageFilter2D::initialise_keymap() +{ + base_type::initialise_keymap(); + this->parser.add_start_key("Gamma Filter Parameters"); + this->parser.add_stop_key("END Gamma Filter Parameters"); +} + +template <> +const char* const GammaImageFilter2D::registered_name = "Gamma"; + +#ifdef _MSC_VER +# pragma warning(disable : 4660) +#endif + +// Explicit instantiation +template class GammaImageFilter2D; + +END_NAMESPACE_STIR diff --git a/src/buildblock/ProjDataInMemory.cxx b/src/buildblock/ProjDataInMemory.cxx index 7f51eb9976..4571e2a8ae 100644 --- a/src/buildblock/ProjDataInMemory.cxx +++ b/src/buildblock/ProjDataInMemory.cxx @@ -356,6 +356,12 @@ ProjDataInMemory::ProjDataInMemory(const ProjDataInMemory& proj_data) std::copy(proj_data.begin_all(), proj_data.end_all(), this->begin_all()); } +shared_ptr +ProjDataInMemory::read_from_file(const std::string& filename) +{ + return std::make_shared(*ProjData::read_from_file(filename)); +} + float ProjDataInMemory::get_bin_value(Bin& bin) { diff --git a/src/buildblock/Scanner.cxx b/src/buildblock/Scanner.cxx index 2f48f14387..4cf34c5c6d 100644 --- a/src/buildblock/Scanner.cxx +++ b/src/buildblock/Scanner.cxx @@ -174,7 +174,7 @@ Scanner::Scanner(Type scanner_type) 7.0F, 6.75F, 3.12932F, - static_cast(15. * _PI / 180), + static_cast(-15. * _PI / 180), 1, 4, 8, @@ -201,7 +201,7 @@ Scanner::Scanner(Type scanner_type) 7.0F, 6.75F, 3.375F, - static_cast(15. * _PI / 180), + static_cast(-15. * _PI / 180), 1, 4, 8, @@ -228,7 +228,7 @@ Scanner::Scanner(Type scanner_type) 7.0F, 6.75F, 3.375F, - static_cast(15. * _PI / 180), + static_cast(-15. * _PI / 180), 3, 4, 8, @@ -255,7 +255,7 @@ Scanner::Scanner(Type scanner_type) 7.0F, 6.25F, 1.650F, - static_cast(13. * _PI / 180), + static_cast(-13. * _PI / 180), 1, 8, 8, @@ -864,9 +864,10 @@ Scanner::Scanner(Type scanner_type) 1, 0.0944F, // energy resolution from Hsu et al. 2017 511.F, - (short int)(0), - (float)(0), // TODO - (float)(0)); + (short int)(2 * 188 + 1), // from RDF file + (float)(13), // from RDF file + (float)(375.4) // Hsu et al. + ); break; case DiscoveryMI4ring: // This is the 4-ring DMI @@ -895,9 +896,10 @@ Scanner::Scanner(Type scanner_type) 1, 0.0944F, // energy resolution from Hsu et al. 2017 511.F, - (short int)(0), - (float)(0), // TODO - (float)(0)); + (short int)(2 * 188 + 1), // from RDF file + (float)(13), // from RDF file + (float)(375.4) // Hsu et al. + ); break; case DiscoveryMI5ring: // This is the 5-ring DMI @@ -926,9 +928,42 @@ Scanner::Scanner(Type scanner_type) 1, 0.0944F, // energy resolution from Hsu et al. 2017 511.F, - (short int)(0), - (float)(0), // TODO - (float)(0)); + (short int)(2 * 188 + 1), // from RDF file + (float)(13), // from RDF file + (float)(375.4) // Hsu et al. + ); + break; + + case DiscoveryMI6ring: // This is the 6-ring DMI + // as above, but one extra block + // Hsu et al. 2017 JNM + // crystal size 3.95 x 5.3 x 25 + set_params(DiscoveryMI6ring, + string_list("GE Discovery MI 6 rings", + "Discovery MI6", + "Discovery MI"), // needs to include last value as used by GE in RDF files + 54, + 415, + 401, // TODO should compute num_arccorrected_bins from effective_FOV/default_bin_size + 2 * 272, + 380.5F - 9.4F, // TODO inner_ring_radius and DOI, currently set such that effective ring-radius is correct + 9.4F, // TODO DOI + 5.52296F, // ring-spacing + 2.206F, // TODO currently using the central bin size default bin size. GE might be using something else + static_cast(-4.399 * _PI / 180), // TODO check sign + 6, + 4, + 9, + 4, + 1, + 1, + 1, + 0.0944F, // energy resolution from Hsu et al. 2017 + 511.F, + (short int)(2 * 188 + 1), // from RDF file + (float)(13), // from RDF file + (float)(375.4) // Hsu et al. + ); break; case HZLR: @@ -2171,6 +2206,7 @@ Scanner::get_scanner_from_name(const string& name) type = static_cast(int_type); } // it's not in the list + warning(std::string("Scanner::get_scanner_from_name: scanner'") + name + "' not found"); return new Scanner(Unknown_scanner); } diff --git a/src/buildblock/WienerArrayFilter2D.cxx b/src/buildblock/WienerArrayFilter2D.cxx new file mode 100644 index 0000000000..8b0e35c730 --- /dev/null +++ b/src/buildblock/WienerArrayFilter2D.cxx @@ -0,0 +1,97 @@ +// +// +/*! + + \file + \ingroup Array + \brief Implementations for class stir::WienerArrayFilter2D + + \author Dimitra Kyriakopoulou + \author Kris Thielemans + +*/ +/* + Copyright (C) 2024, University College London + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ +#include "stir/WienerArrayFilter2D.h" + +#include + +using std::nth_element; + +START_NAMESPACE_STIR + +template +WienerArrayFilter2D::WienerArrayFilter2D() +{} + +template +void +WienerArrayFilter2D::do_it(Array<3, elemT>& out_array, const Array<3, elemT>& in_array) const +{ + assert(out_array.get_index_range() == in_array.get_index_range()); + + const int sx = out_array[0][0].get_length(); + const int sy = out_array[0].get_length(); + const int sa = out_array.get_length(); + const int min_x = in_array[0][0].get_min_index(); + const int min_y = in_array[0].get_min_index(); + const int ws = 9; + + for (int ia = 0; ia < sa; ia++) + { + std::vector> localMean(sx, std::vector(sy, 0.0f)); + std::vector> localVar(sx, std::vector(sy, 0.0f)); + + float noise = 0.; + + for (int i = 1; i < sx - 1; i++) + { + for (int j = 1; j < sy - 1; j++) + { + localMean[i][j] = 0; + localVar[i][j] = 0; + + for (int k = -1; k <= 1; k++) + for (int l = -1; l <= 1; l++) + localMean[i][j] += in_array[ia][min_x + i + k][min_y + j + l]; + localMean[i][j] /= ws; + + for (int k = -1; k <= 1; k++) + for (int l = -1; l <= 1; l++) + localVar[i][j] += in_array[ia][min_x + i + k][min_y + j + l] * in_array[ia][min_x + i + k][min_y + j + l]; + localVar[i][j] = localVar[i][j] / ws - localMean[i][j] * localMean[i][j]; + + noise += localVar[i][j]; + } + } + noise /= (sx * sy); + + for (int i = 1; i < sx - 1; i++) + { + for (int j = 1; j < sy - 1; j++) + { + out_array[ia][min_x + i][min_y + j] = (in_array[ia][min_x + i][min_y + j] - localMean[i][j]) + / std::max(localVar[i][j], noise) * std::max(localVar[i][j] - noise, 0.f) + + localMean[i][j]; + } + } + } +} + +template +bool +WienerArrayFilter2D::is_trivial() const +{ + return false; +} + +// instantiation +template class WienerArrayFilter2D; + +END_NAMESPACE_STIR diff --git a/src/buildblock/WienerImageFilter2D.cxx b/src/buildblock/WienerImageFilter2D.cxx new file mode 100644 index 0000000000..96fcaca2da --- /dev/null +++ b/src/buildblock/WienerImageFilter2D.cxx @@ -0,0 +1,94 @@ +// +// +/*! + \file + \ingroup ImageProcessor + \brief Implementations for class stir::WienerImageFilter2D + + \author Dimitra Kyriakopoulou + \author Kris Thielemans + +*/ +/* + Copyright (C) 2024, University College London + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +#include "stir/WienerImageFilter2D.h" +#include "stir/CartesianCoordinate2D.h" +#include "stir/DiscretisedDensity.h" + +START_NAMESPACE_STIR + +template +WienerImageFilter2D::WienerImageFilter2D() +{ + std::cout << "Wiener filter start" << std::endl; + + set_defaults(); +} + +template +Succeeded +WienerImageFilter2D::virtual_set_up(const DiscretisedDensity<3, elemT>& density) +{ + + /* if (consistency_check(density) == Succeeded::no) + return Succeeded::no;*/ + wiener_filter = WienerArrayFilter2D(); + + return Succeeded::yes; +} + +template +void +WienerImageFilter2D::virtual_apply(DiscretisedDensity<3, elemT>& density) const +{ + // assert(consistency_check(density) == Succeeded::yes); + wiener_filter(density); +} + +template +void +WienerImageFilter2D::virtual_apply(DiscretisedDensity<3, elemT>& out_density, + const DiscretisedDensity<3, elemT>& in_density) const +{ + // assert(consistency_check(in_density) == Succeeded::yes); + wiener_filter(out_density, in_density); +} + +template +void +WienerImageFilter2D::set_defaults() +{ + base_type::set_defaults(); +} + +template +void +WienerImageFilter2D::initialise_keymap() +{ + base_type::initialise_keymap(); + this->parser.add_start_key("Wiener Filter Parameters"); + this->parser.add_stop_key("END Wiener Filter Parameters"); +} + +template <> +const char* const WienerImageFilter2D::registered_name = "Wiener"; + +#ifdef _MSC_VER +// prevent warning message on reinstantiation, +// note that we get a linking error if we don't have the explicit instantiation below +# pragma warning(disable : 4660) +#endif + +// Register this class in the ImageProcessor registry +// static WienerImageFilter2D::RegisterIt dummyWiener; + +template class WienerImageFilter2D; + +END_NAMESPACE_STIR diff --git a/src/buildblock/buildblock_registries.cxx b/src/buildblock/buildblock_registries.cxx index 577f0b1a42..1474f6cc66 100644 --- a/src/buildblock/buildblock_registries.cxx +++ b/src/buildblock/buildblock_registries.cxx @@ -2,7 +2,7 @@ // /* Copyright (C) 2000- 2007, Hammersmith Imanet Ltd - Copyright (C) 2019, 2020, UCL + Copyright (C) 2019, 2020, 2024 UCL This file is part of STIR. SPDX-License-Identifier: Apache-2.0 @@ -17,12 +17,14 @@ \brief File that registers all stir::RegisterObject children in buildblock \author Kris Thielemans - + \author Dimitra Kyriakopoulou */ #include "stir/SeparableCartesianMetzImageFilter.h" #include "stir/SeparableGaussianImageFilter.h" #include "stir/MedianImageFilter3D.h" +#include "stir/WienerImageFilter2D.h" +#include "stir/GammaImageFilter2D.h" #include "stir/MinimalImageFilter3D.h" #include "stir/ChainedDataProcessor.h" #include "stir/ThresholdMinToSmallPositiveValueDataProcessor.h" @@ -35,6 +37,8 @@ START_NAMESPACE_STIR static MedianImageFilter3D::RegisterIt dummy; +static WienerImageFilter2D::RegisterIt dummyWiener; +static GammaImageFilter2D::RegisterIt dummyGamma; static MinimalImageFilter3D::RegisterIt dummy1; static SeparableCartesianMetzImageFilter::RegisterIt dummy2; static SeparableGaussianImageFilter::RegisterIt dummySGF; diff --git a/src/buildblock/extend_projdata.cxx b/src/buildblock/extend_projdata.cxx index 1ecd11c158..59fa14288c 100644 --- a/src/buildblock/extend_projdata.cxx +++ b/src/buildblock/extend_projdata.cxx @@ -93,7 +93,7 @@ extend_segment(const SegmentBySinogram& segment, if (extend_without_wrapping) { out[axial_pos][min_dim[2] + view_edge] = out[axial_pos][min_dim[2] + view_extension]; - out[axial_pos][max_dim[2] - view_extension] = out[axial_pos][max_dim[2] - view_extension]; + out[axial_pos][max_dim[2] - view_edge] = out[axial_pos][max_dim[2] - view_extension]; } else if (flip_views) { diff --git a/src/buildblock/interpolate_projdata.cxx b/src/buildblock/interpolate_projdata.cxx index 060e8aa35d..8e0eee6a87 100644 --- a/src/buildblock/interpolate_projdata.cxx +++ b/src/buildblock/interpolate_projdata.cxx @@ -23,6 +23,7 @@ //#include "stir/display.h" #include "stir/ProjDataInfo.h" #include "stir/ProjDataInfoCylindricalNoArcCorr.h" +#include "stir/ProjDataInfoGenericNoArcCorr.h" #include "stir/IndexRange.h" #include "stir/BasicCoordinate.h" #include "stir/Sinogram.h" @@ -191,6 +192,10 @@ interpolate_projdata(ProjData& proj_data_out, error("interpolate_projdata needs both projection to be of a scanner with the same ring radius"); } + // handle BlocksOnCylindrical interpolation manually + if (proj_data_in_info.get_scanner_sptr()->get_scanner_geometry() != "Cylindrical") + return interpolate_blocks_on_cylindrical_projdata(proj_data_out, proj_data_in, remove_interleaving); + // initialise interpolator BSpline::BSplinesRegularGrid<3, float, float> proj_data_interpolator(these_types); for (int k = proj_data_out_info.get_min_tof_pos_num(); k <= proj_data_out_info.get_max_tof_pos_num(); ++k) @@ -200,96 +205,319 @@ interpolate_projdata(ProjData& proj_data_out, proj_data_in.get_segment_by_sinogram(0, k)) : proj_data_in.get_segment_by_sinogram(0, k); + // for Cylindrical, spacing is regular in all directions, which makes mapping trivial std::function(const BasicCoordinate<3, int>&)> index_converter; - if (proj_data_in_info.get_scanner_sptr()->get_scanner_geometry() == "Cylindrical") - { // for Cylindrical, spacing is regular in all directions, which makes mapping trivial - // especially in view direction, extending by 5 leads to much smaller artifacts - proj_data_interpolator.set_coef(extend_segment(segment, 5, 5, 5)); - - BasicCoordinate<3, double> offset, step; - // out_index * step + offset = in_index - const float in_sampling_m = proj_data_in_info.get_sampling_in_m(Bin(0, 0, 0, 0)); - const float out_sampling_m = proj_data_out_info.get_sampling_in_m(Bin(0, 0, 0, 0)); - // offset in 'in' index units - offset[1] = (proj_data_out_info.get_m(Bin(0, 0, 0, 0)) - proj_data_in_info.get_m(Bin(0, 0, 0, 0))) / in_sampling_m; - step[1] = out_sampling_m / in_sampling_m; - - const float in_sampling_phi = (proj_data_in_info.get_phi(Bin(0, 1, 0, 0)) - proj_data_in_info.get_phi(Bin(0, 0, 0, 0))) - / (remove_interleaving ? 2 : 1); - const float out_sampling_phi - = proj_data_out_info.get_phi(Bin(0, 1, 0, 0)) - proj_data_out_info.get_phi(Bin(0, 0, 0, 0)); - offset[2] - = (proj_data_out_info.get_phi(Bin(0, 0, 0, 0)) - proj_data_in_info.get_phi(Bin(0, 0, 0, 0))) / in_sampling_phi; - step[2] = out_sampling_phi / in_sampling_phi; - - const float in_sampling_s = proj_data_in_info.get_sampling_in_s(Bin(0, 0, 0, 0)); - const float out_sampling_s = proj_data_out_info.get_sampling_in_s(Bin(0, 0, 0, 0)); - offset[3] = (proj_data_out_info.get_s(Bin(0, 0, 0, 0)) - proj_data_in_info.get_s(Bin(0, 0, 0, 0))) / in_sampling_s; - step[3] = out_sampling_s / in_sampling_s; - - // define a function to translate indices in the output proj data to indices in input proj data - index_converter - = [&proj_data_out_info, offset, step](const BasicCoordinate<3, int>& index_out) -> BasicCoordinate<3, double> { - // translate to indices in input proj data - BasicCoordinate<3, double> index_in; - for (auto dim = 1; dim <= 3; dim++) - index_in[dim] = index_out[dim] * step[dim] + offset[dim]; - - return index_in; - }; - } - else - { // for BlocksOnCylindrical, views and tangential positions are not subsampled and can be mapped 1:1 - if (proj_data_in_info.get_num_tangential_poss() != proj_data_out_info.get_num_tangential_poss()) - { - error("Interpolation of BlocksOnCylindrical scanners assumes that number of tangential positions " - "is the same in the downsampled scanner."); - } - if (proj_data_in_info.get_num_views() != proj_data_out_info.get_num_views()) - { - error("Interpolation of BlocksOnCylindrical scanners assumes that number of views " - "is the same in the downsampled scanner."); - } + // especially in view direction, extending by 5 leads to much smaller artifacts + proj_data_interpolator.set_coef(extend_segment(segment, 5, 5, 5)); + + BasicCoordinate<3, double> offset, step; + // out_index * step + offset = in_index + const float in_sampling_m = proj_data_in_info.get_sampling_in_m(Bin(0, 0, 0, 0)); + const float out_sampling_m = proj_data_out_info.get_sampling_in_m(Bin(0, 0, 0, 0)); + // offset in 'in' index units + offset[1] = (proj_data_out_info.get_m(Bin(0, 0, 0, 0)) - proj_data_in_info.get_m(Bin(0, 0, 0, 0))) / in_sampling_m; + step[1] = out_sampling_m / in_sampling_m; + + const float in_sampling_phi = (proj_data_in_info.get_phi(Bin(0, 1, 0, 0)) - proj_data_in_info.get_phi(Bin(0, 0, 0, 0))) + / (remove_interleaving ? 2 : 1); + const float out_sampling_phi = proj_data_out_info.get_phi(Bin(0, 1, 0, 0)) - proj_data_out_info.get_phi(Bin(0, 0, 0, 0)); + offset[2] = (proj_data_out_info.get_phi(Bin(0, 0, 0, 0)) - proj_data_in_info.get_phi(Bin(0, 0, 0, 0))) / in_sampling_phi; + step[2] = out_sampling_phi / in_sampling_phi; + + const float in_sampling_s = proj_data_in_info.get_sampling_in_s(Bin(0, 0, 0, 0)); + const float out_sampling_s = proj_data_out_info.get_sampling_in_s(Bin(0, 0, 0, 0)); + offset[3] = (proj_data_out_info.get_s(Bin(0, 0, 0, 0)) - proj_data_in_info.get_s(Bin(0, 0, 0, 0))) / in_sampling_s; + step[3] = out_sampling_s / in_sampling_s; + + // define a function to translate indices in the output proj data to indices in input proj data + index_converter + = [&proj_data_out_info, offset, step](const BasicCoordinate<3, int>& index_out) -> BasicCoordinate<3, double> { + // translate to indices in input proj data + BasicCoordinate<3, double> index_in; + for (auto dim = 1; dim <= 3; dim++) + index_in[dim] = index_out[dim] * step[dim] + offset[dim]; + + return index_in; + }; - // only extending in axial direction - an extension of 2 was found to be sufficient - proj_data_interpolator.set_coef(extend_segment(segment, 0, 2, 0)); + SegmentBySinogram sino_3D_out = proj_data_out.get_empty_segment_by_sinogram(0, false, k); + sample_function_using_index_converter(sino_3D_out, proj_data_interpolator, index_converter); - auto m_offset = proj_data_in_info.get_m(Bin(0, 0, 0, 0)); - auto m_sampling = proj_data_in_info.get_sampling_in_m(Bin(0, 0, 0, 0)); + if (proj_data_out.set_segment(sino_3D_out) == Succeeded::no) + return Succeeded::no; + } + return Succeeded::yes; +} - // confirm that proj_data_in has equidistant sampling in m - for (auto axial_pos = proj_data_in_info.get_min_axial_pos_num(0); - axial_pos <= proj_data_in_info.get_max_axial_pos_num(0); - axial_pos++) - { - if (abs(m_sampling - proj_data_in_info.get_sampling_in_m(Bin(0, 0, axial_pos, 0))) > 1E-4) - error("input projdata to interpolate_projdata are not equidistantly sampled in m."); - } +//! This function interpolates BlocksOnCylindrical proj data taking bucket intersections and gaps into account. +/*! The above proj data interpolation function works well for cylindrical scanners where there are no large discontinuities. +For BlocksOnCylindrical scanners, the sudden change of angle from one bucket to the next and the gaps between blocks and buckets +lead to interpolation artifacts if not taken into account. Therefore, this function does the interpolation in physical space +rather than directly on the proj data bin values. For each bin in proj_data_out (the full size proj data), we find the four +closest LORs in the downsampled proj_data_in. These are then weighted using bilinear interpolation based on the crystal positions +of the two endpoints of the LOR. +*/ +Succeeded +interpolate_blocks_on_cylindrical_projdata(ProjData& proj_data_out, const ProjData& proj_data_in, bool remove_interleaving) +{ + const ProjDataInfo& proj_data_in_info = *proj_data_in.get_proj_data_info_sptr(); + const ProjDataInfo& proj_data_out_info = *proj_data_out.get_proj_data_info_sptr(); - // define a function to translate indices in the output proj data to indices in input proj data - index_converter = [&proj_data_out_info, m_offset, m_sampling]( - const BasicCoordinate<3, int>& index_out) -> BasicCoordinate<3, double> { - // translate index on output to coordinate - auto bin - = Bin(0 /* segment */, index_out[2] /* view */, index_out[1] /* axial pos */, index_out[3] /* tangential pos */); - auto out_m = proj_data_out_info.get_m(bin); - - // translate to indices in input proj data - BasicCoordinate<3, double> index_in; - index_in[1] = (out_m - m_offset) / m_sampling; - index_in[2] = index_out[2]; - index_in[3] = index_out[3]; - - return index_in; - }; - } + auto m_offset = proj_data_in_info.get_m(Bin(0, 0, 0, 0)); + auto m_sampling = proj_data_in_info.get_sampling_in_m(Bin(0, 0, 0, 0)); + for (int k = proj_data_out_info.get_min_tof_pos_num(); k <= proj_data_out_info.get_max_tof_pos_num(); ++k) + { + SegmentBySinogram segment + = remove_interleaving ? make_non_interleaved_segment(*(make_non_interleaved_proj_data_info(proj_data_in_info)), + proj_data_in.get_segment_by_sinogram(0, k)) + : proj_data_in.get_segment_by_sinogram(0, k); SegmentBySinogram sino_3D_out = proj_data_out.get_empty_segment_by_sinogram(0, false, k); - sample_function_using_index_converter(sino_3D_out, proj_data_interpolator, index_converter); + BasicCoordinate<3, int> min_out, max_out; + IndexRange<3> out_range = sino_3D_out.get_index_range(); + if (!out_range.get_regular_range(min_out, max_out)) + warning("Output must be regular range!"); + + // confirm that proj_data_in has equidistant sampling in m + for (auto axial_pos = proj_data_in_info.get_min_axial_pos_num(0); axial_pos <= proj_data_in_info.get_max_axial_pos_num(0); + axial_pos++) + { + if (abs(m_sampling - proj_data_in_info.get_sampling_in_m(Bin(0, 0, axial_pos, 0))) > 1E-4) + error("input projdata to interpolate_projdata are not equidistantly sampled in m."); + } + + BasicCoordinate<3, int> index_out; + for (index_out[1] = min_out[1]; index_out[1] <= max_out[1]; ++index_out[1]) + { + for (index_out[2] = min_out[2]; index_out[2] <= max_out[2]; ++index_out[2]) + { + for (index_out[3] = min_out[3]; index_out[3] <= max_out[3]; ++index_out[3]) + { + // translate index on output to coordinate + auto bin + = Bin(0 /* segment */, index_out[2] /* view */, + index_out[1] /* axial pos */, index_out[3] /* tangential pos + */); + auto out_m = proj_data_out_info.get_m(bin); + double axial_idx = (out_m - m_offset) / m_sampling; + int axial_floor = std::floor(axial_idx); + int axial_ceil = std::ceil(axial_idx); + if (axial_floor == axial_ceil) + { + if (axial_floor == 0) + axial_ceil++; + else + axial_floor--; + } + + // find the two crystals for this bin and on which buckets (modules) they are + int det1_num_out, det2_num_out; + const auto proj_data_out_info_ptr = dynamic_cast(&proj_data_out_info); + proj_data_out_info_ptr->get_det_num_pair_for_view_tangential_pos_num(det1_num_out, + det2_num_out, + index_out[2], /*view*/ + index_out[3] /* tangential pos */); + const int dets_per_module_out = proj_data_out_info.get_scanner_sptr()->get_num_transaxial_crystals_per_bucket(); + const int det1_module = det1_num_out / dets_per_module_out; + const int det2_module = det2_num_out / dets_per_module_out; + + // translate the crystal position on each module from the full size scanner to the downsampled scanner + const auto proj_data_in_info_ptr = dynamic_cast(&proj_data_in_info); + const int dets_per_module_in + = proj_data_in_info_ptr->get_scanner_sptr()->get_num_transaxial_crystals_per_bucket(); + const int crystal1_out_module_idx = det1_num_out % dets_per_module_out; + const double crystal1_out_module_pos + = std::floor(static_cast(crystal1_out_module_idx) + / proj_data_out_info_ptr->get_scanner_sptr()->get_num_transaxial_crystals_per_block()) + * proj_data_out_info_ptr->get_scanner_sptr()->get_transaxial_block_spacing() + + static_cast( + crystal1_out_module_idx + % proj_data_out_info_ptr->get_scanner_sptr()->get_num_transaxial_crystals_per_block()) + * proj_data_out_info_ptr->get_scanner_sptr()->get_transaxial_crystal_spacing(); + const double crystal1_num_in + = det1_module * dets_per_module_in + + crystal1_out_module_pos / proj_data_in_info.get_scanner_sptr()->get_transaxial_crystal_spacing(); + const int crystal2_out_module_idx = det2_num_out % dets_per_module_out; + const double crystal2_out_module_pos + = std::floor(static_cast(crystal2_out_module_idx) + / proj_data_out_info_ptr->get_scanner_sptr()->get_num_transaxial_crystals_per_block()) + * proj_data_out_info_ptr->get_scanner_sptr()->get_transaxial_block_spacing() + + static_cast( + crystal2_out_module_idx + % proj_data_out_info_ptr->get_scanner_sptr()->get_num_transaxial_crystals_per_block()) + * proj_data_out_info_ptr->get_scanner_sptr()->get_transaxial_crystal_spacing(); + const double crystal2_num_in + = det2_module * dets_per_module_in + + crystal2_out_module_pos / proj_data_in_info.get_scanner_sptr()->get_transaxial_crystal_spacing(); + const auto crystal1_num_in_floor + = std::max(static_cast(std::floor(crystal1_num_in)), det1_module * dets_per_module_in); + const auto crystal1_num_in_ceil + = std::min(static_cast(std::ceil(crystal1_num_in)), (det1_module + 1) * dets_per_module_in - 1); + const auto crystal2_num_in_floor + = std::max(static_cast(std::floor(crystal2_num_in)), det2_module * dets_per_module_in); + const auto crystal2_num_in_ceil + = std::min(static_cast(std::ceil(crystal2_num_in)), (det2_module + 1) * dets_per_module_in - 1); + + // translate to indices in input proj data + BasicCoordinate<3, int> index_in; + + // in this case we can skip parts of the interpolation + if (crystal1_num_in_floor == crystal1_num_in_ceil) + { + if (crystal2_num_in_floor == crystal2_num_in_ceil) + { + int one_view, one_tang; + proj_data_in_info_ptr->get_view_tangential_pos_num_for_det_num_pair( + one_view, one_tang, crystal1_num_in_floor, crystal2_num_in_floor); + one_tang = std::min(std::max(proj_data_in_info_ptr->get_min_tangential_pos_num(), one_tang), + proj_data_in_info_ptr->get_max_tangential_pos_num()); + + index_in[1] = std::max(axial_floor, proj_data_in_info_ptr->get_min_axial_pos_num(0)); + index_in[2] = one_view; + index_in[3] = one_tang; + sino_3D_out[index_out] += segment[index_in] * (axial_ceil - axial_idx); + + index_in[1] = std::min(axial_ceil, proj_data_in_info_ptr->get_max_axial_pos_num(0)); + sino_3D_out[index_out] += segment[index_in] * (axial_idx - axial_floor); + } + else + { + int ff_view, fc_view; + int ff_tang, fc_tang; + proj_data_in_info_ptr->get_view_tangential_pos_num_for_det_num_pair( + ff_view, ff_tang, crystal1_num_in_floor, crystal2_num_in_floor); + proj_data_in_info_ptr->get_view_tangential_pos_num_for_det_num_pair( + fc_view, fc_tang, crystal1_num_in_floor, crystal2_num_in_ceil); + + ff_tang = std::min(std::max(proj_data_in_info_ptr->get_min_tangential_pos_num(), ff_tang), + proj_data_in_info_ptr->get_max_tangential_pos_num()); + fc_tang = std::min(std::max(proj_data_in_info_ptr->get_min_tangential_pos_num(), fc_tang), + proj_data_in_info_ptr->get_max_tangential_pos_num()); + + index_in[1] = std::max(axial_floor, proj_data_in_info_ptr->get_min_axial_pos_num(0)); + index_in[2] = ff_view; + index_in[3] = ff_tang; + sino_3D_out[index_out] + += segment[index_in] * (axial_ceil - axial_idx) * (crystal2_num_in_ceil - crystal2_num_in); + index_in[2] = fc_view; + index_in[3] = fc_tang; + sino_3D_out[index_out] + += segment[index_in] * (axial_ceil - axial_idx) * (crystal2_num_in - crystal2_num_in_floor); + index_in[1] = std::min(axial_ceil, proj_data_in_info_ptr->get_max_axial_pos_num(0)); + sino_3D_out[index_out] + += segment[index_in] * (axial_idx - axial_floor) * (crystal2_num_in - crystal2_num_in_floor); + index_in[2] = ff_view; + index_in[3] = ff_tang; + sino_3D_out[index_out] + += segment[index_in] * (axial_idx - axial_floor) * (crystal2_num_in_ceil - crystal2_num_in); + } + } + else if (crystal2_num_in_floor == crystal2_num_in_ceil) + { + int ff_view, cf_view; + int ff_tang, cf_tang; + proj_data_in_info_ptr->get_view_tangential_pos_num_for_det_num_pair( + ff_view, ff_tang, crystal1_num_in_floor, crystal2_num_in_floor); + proj_data_in_info_ptr->get_view_tangential_pos_num_for_det_num_pair( + cf_view, cf_tang, crystal1_num_in_ceil, crystal2_num_in_floor); + + ff_tang = std::min(std::max(proj_data_in_info_ptr->get_min_tangential_pos_num(), ff_tang), + proj_data_in_info_ptr->get_max_tangential_pos_num()); + cf_tang = std::min(std::max(proj_data_in_info_ptr->get_min_tangential_pos_num(), cf_tang), + proj_data_in_info_ptr->get_max_tangential_pos_num()); + + index_in[1] = std::max(axial_floor, proj_data_in_info_ptr->get_min_axial_pos_num(0)); + index_in[2] = ff_view; + index_in[3] = ff_tang; + sino_3D_out[index_out] + += segment[index_in] * (axial_ceil - axial_idx) * (crystal1_num_in_ceil - crystal1_num_in); + index_in[2] = cf_view; + index_in[3] = cf_tang; + sino_3D_out[index_out] + += segment[index_in] * (axial_ceil - axial_idx) * (crystal1_num_in - crystal1_num_in_floor); + index_in[1] = std::min(axial_ceil, proj_data_in_info_ptr->get_max_axial_pos_num(0)); + sino_3D_out[index_out] + += segment[index_in] * (axial_idx - axial_floor) * (crystal1_num_in - crystal1_num_in_floor); + index_in[2] = ff_view; + index_in[3] = ff_tang; + sino_3D_out[index_out] + += segment[index_in] * (axial_idx - axial_floor) * (crystal1_num_in_ceil - crystal1_num_in); + } + else // in this case we need to do a trilinear interpolation + { + int ff_view, fc_view, cf_view, cc_view; + int ff_tang, fc_tang, cf_tang, cc_tang; + proj_data_in_info_ptr->get_view_tangential_pos_num_for_det_num_pair( + ff_view, ff_tang, crystal1_num_in_floor, crystal2_num_in_floor); + proj_data_in_info_ptr->get_view_tangential_pos_num_for_det_num_pair( + fc_view, fc_tang, crystal1_num_in_floor, crystal2_num_in_ceil); + proj_data_in_info_ptr->get_view_tangential_pos_num_for_det_num_pair( + cf_view, cf_tang, crystal1_num_in_ceil, crystal2_num_in_floor); + proj_data_in_info_ptr->get_view_tangential_pos_num_for_det_num_pair( + cc_view, cc_tang, crystal1_num_in_ceil, crystal2_num_in_ceil); + + // TODO: why can we get positions out that are not even in the proj data?! + ff_tang = std::min(std::max(proj_data_in_info_ptr->get_min_tangential_pos_num(), ff_tang), + proj_data_in_info_ptr->get_max_tangential_pos_num()); + fc_tang = std::min(std::max(proj_data_in_info_ptr->get_min_tangential_pos_num(), fc_tang), + proj_data_in_info_ptr->get_max_tangential_pos_num()); + cf_tang = std::min(std::max(proj_data_in_info_ptr->get_min_tangential_pos_num(), cf_tang), + proj_data_in_info_ptr->get_max_tangential_pos_num()); + cc_tang = std::min(std::max(proj_data_in_info_ptr->get_min_tangential_pos_num(), cc_tang), + proj_data_in_info_ptr->get_max_tangential_pos_num()); + + index_in[1] = std::max(axial_floor, proj_data_in_info_ptr->get_min_axial_pos_num(0)); + index_in[2] = ff_view; + index_in[3] = ff_tang; + sino_3D_out[index_out] += segment[index_in] * (axial_ceil - axial_idx) + * (crystal1_num_in_ceil - crystal1_num_in) + * (crystal2_num_in_ceil - crystal2_num_in); + index_in[2] = fc_view; + index_in[3] = fc_tang; + sino_3D_out[index_out] += segment[index_in] * (axial_ceil - axial_idx) + * (crystal1_num_in_ceil - crystal1_num_in) + * (crystal2_num_in - crystal2_num_in_floor); + index_in[2] = cc_view; + index_in[3] = cc_tang; + sino_3D_out[index_out] += segment[index_in] * (axial_ceil - axial_idx) + * (crystal1_num_in - crystal1_num_in_floor) + * (crystal2_num_in - crystal2_num_in_floor); + index_in[2] = cf_view; + index_in[3] = cf_tang; + sino_3D_out[index_out] += segment[index_in] * (axial_ceil - axial_idx) + * (crystal1_num_in - crystal1_num_in_floor) + * (crystal2_num_in_ceil - crystal2_num_in); + + index_in[1] = std::min(axial_ceil, proj_data_in_info_ptr->get_max_axial_pos_num(0)); + index_in[2] = ff_view; + index_in[3] = ff_tang; + sino_3D_out[index_out] += segment[index_in] * (axial_idx - axial_floor) + * (crystal1_num_in_ceil - crystal1_num_in) + * (crystal2_num_in_ceil - crystal2_num_in); + index_in[2] = fc_view; + index_in[3] = fc_tang; + sino_3D_out[index_out] += segment[index_in] * (axial_idx - axial_floor) + * (crystal1_num_in_ceil - crystal1_num_in) + * (crystal2_num_in - crystal2_num_in_floor); + index_in[2] = cc_view; + index_in[3] = cc_tang; + sino_3D_out[index_out] += segment[index_in] * (axial_idx - axial_floor) + * (crystal1_num_in - crystal1_num_in_floor) + * (crystal2_num_in - crystal2_num_in_floor); + index_in[2] = cf_view; + index_in[3] = cf_tang; + sino_3D_out[index_out] += segment[index_in] * (axial_idx - axial_floor) + * (crystal1_num_in - crystal1_num_in_floor) + * (crystal2_num_in_ceil - crystal2_num_in); + } + } + } + } if (proj_data_out.set_segment(sino_3D_out) == Succeeded::no) return Succeeded::no; } + return Succeeded::yes; } diff --git a/src/cmake/FindCERN_ROOT.cmake b/src/cmake/FindCERN_ROOT.cmake index 4879779640..9aad3924f4 100644 --- a/src/cmake/FindCERN_ROOT.cmake +++ b/src/cmake/FindCERN_ROOT.cmake @@ -14,8 +14,9 @@ ## FINDING ROOT # # Attempts to `find_package(ROOT)`, If that fails, use root-config. -# The primary method for ROOT being found is to use the `find_package(ROOT ${CERN_ROOT_FIND_VERSION} QUIET)` call. -# This process utilizes the `ROOT_DIR` variable to find the relevant CMake files. +# The primary method for ROOT being found is to use the `find_package(ROOT ${CERN_ROOT_FIND_VERSION} CONFIG)` call. +# This process utilizes the `ROOT_DIR` variable to find the relevant ROOTConfig*.cmake files (which are +# part of the ROOT distribution). # There are two methods by which this variable can be set: # 1. Set the `ROOT_DIR` CMake argument. Point to `/cmake` directory. # 2. Use the `ROOTSYS` CMake or environment variable. If `ROOT_DIR` is not provided, we will determine set `ROOT_DIR` to `${ROOTSYS}/cmake`. @@ -53,7 +54,7 @@ if (DEFINED ROOT_DIR) endif() endif() -find_package(ROOT ${CERN_ROOT_FIND_VERSION} QUIET) +find_package(ROOT ${CERN_ROOT_FIND_VERSION} CONFIG) if (ROOT_FOUND) if (CERN_ROOT_DEBUG) message(STATUS "Found ROOTConfig.cmake, so translating to old CERN_ROOT variable names") diff --git a/src/cmake/STIRConfig.cmake.in b/src/cmake/STIRConfig.cmake.in index 52fd60b592..0b1851b06d 100644 --- a/src/cmake/STIRConfig.cmake.in +++ b/src/cmake/STIRConfig.cmake.in @@ -63,15 +63,24 @@ else() endif() ## find external packages +# Note that we generally need find_package even for those that are only +# PRIVATE dependencies (such as ITK, parallelproj etc). This is because +# otherwise their targets (such as parallelproj::parallelproj_c) don't exist, +# and the user will get CMake errors. +# See https://github.com/UCL/STIR/issues/1535 + # we use a trick by Matt McCormick (kitware) # to set ITK_DIR etc first before calling find_package # to make sure we pick the same version of the external library + if (STIR_FIND_QUIETLY) SET(STIR_FIND_TYPE "QUIET") else() SET(STIR_FIND_TYPE "REQUIRED") endif() +find_package(Boost @Boost_VERSION_STRING@ REQUIRED) + if (@ITK_FOUND@) message(STATUS "ITK support in STIR enabled.") set(ITK_DIR "@ITK_DIR@") @@ -79,7 +88,6 @@ if (@ITK_FOUND@) if(NOT ITK_FOUND) SET(STIR_FOUND OFF) endif() - include(${ITK_USE_FILE}) set(STIR_BUILT_WITH_ITK TRUE) endif() diff --git a/src/cmake/stir_dirs.cmake b/src/cmake/stir_dirs.cmake index c2526f3d9e..f875e4ecea 100644 --- a/src/cmake/stir_dirs.cmake +++ b/src/cmake/stir_dirs.cmake @@ -41,7 +41,7 @@ ${PROJECT_SOURCE_DIR}/src/scatter_buildblock/scatter_registries.cxx ) # need to list IO first such that its dependencies on other libraries are resolved -SET( STIR_LIBRARIES IO analytic_FBP3DRP analytic_FBP2D iterative_OSMAPOSL iterative_KOSMAPOSL +SET( STIR_LIBRARIES IO analytic_FBP3DRP analytic_FBP2D analytic_SRT2D analytic_SRT2DSPECT iterative_OSMAPOSL iterative_KOSMAPOSL iterative_OSSPS scatter_buildblock modelling_buildblock listmode_buildblock recon_buildblock display data_buildblock numerics_buildblock buildblock @@ -71,6 +71,8 @@ SET( STIR_DIRS modelling_utilities listmode_utilities analytic/FBP2D + analytic/SRT2D + analytic/SRT2DSPECT analytic/FBP3DRP iterative/OSMAPOSL iterative/KOSMAPOSL diff --git a/src/include/stir/DetectorCoordinateMap.h b/src/include/stir/DetectorCoordinateMap.h index 5e427a40cd..1278abc1c8 100644 --- a/src/include/stir/DetectorCoordinateMap.h +++ b/src/include/stir/DetectorCoordinateMap.h @@ -123,6 +123,8 @@ class DetectorCoordinateMap explicit DetectorCoordinateMap(double sigma = 0.0) : sigma(sigma) { + if (sigma == 0.0) + return; # ifdef STIR_OPENMP generators.resize(omp_get_max_threads()); distributions.resize(omp_get_max_threads()); diff --git a/src/include/stir/GammaArrayFilter2D.h b/src/include/stir/GammaArrayFilter2D.h new file mode 100644 index 0000000000..1201f4847a --- /dev/null +++ b/src/include/stir/GammaArrayFilter2D.h @@ -0,0 +1,90 @@ +// +// +/*! + \file + \ingroup Array + \brief Declaration of class stir::GammaArrayFilter2D + + \author Dimitra Kyriakopoulou + \author Kris Thielemans +*/ +/* + Copyright (C) 2024, University College London + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +#ifndef __stir_GammaArrayFilter2D_H__ +#define __stir_GammaArrayFilter2D_H__ + +#include "stir/ArrayFunctionObject_2ArgumentImplementation.h" + +START_NAMESPACE_STIR + +/*! + \ingroup Array + \brief Gamma correction filter for 2D slices in a 3D volume. + + This filter enhances image contrast by adjusting pixel intensities using gamma correction. + The algorithm operates on each 2D slice (axial direction) independently and involves: + + 1. **Normalization**: Each pixel value in the slice is normalized to [0, 1]: + \f[ + \text{normalized}(i, j) = \frac{\text{input}(i, j) - \text{min\_val}}{\text{max\_val} - \text{min\_val}} + \f] + where \c min_val and \c max_val are the minimum and maximum pixel values in the slice. + + 2. **Average Pixel Value Calculation**: The filter computes the average pixel value across all pixels in the normalized slice + that have absolute values greater than 0.1. This average value is used to determine the gamma exponent. + + 3. **Gamma Calculation**: Determines the gamma exponent using: + \f[ + \gamma = \frac{\log(0.25)}{\log(\text{averagePixelValue})} + \f] + where \c 0.25 is the target average intensity level for contrast adjustment. + + 4. **Gamma Correction**: Adjusts pixel values using: + \f[ + \text{corrected}(i, j) = \text{normalized}(i, j)^{\gamma} + \f] + + 5. **Rescaling**: Converts the corrected values back to their original range: + \f[ + \text{output}(i, j) = \text{corrected}(i, j) \times (\text{max\_val} - \text{min\_val}) + \text{min\_val} + \f] + + ### Edge Handling + - The filter processes each 2D slice independently and does not apply padding. Therefore, the edges of the slices are fully + processed without boundary effects. + + ### Gamma Filter Configuration + This filter is fully automated and does not require any parameters. + To enable it in the reconstruction process, include the following in the parameter file: + \code + post-filter type := Gamma + Gamma Filter Parameters := + End Gamma Filter Parameters := + \endcode + + \note This filter operates on each axial slice independently, and does not take into account neighboring slices. It is + effectively a 2D filter applied slice-by-slice on a 3D volume. + + */ + +template +class GammaArrayFilter2D : public ArrayFunctionObject_2ArgumentImplementation<3, elemT> +{ +public: + explicit GammaArrayFilter2D(); + bool is_trivial() const override; + +private: + void do_it(Array<3, elemT>& out_array, const Array<3, elemT>& in_array) const override; +}; + +END_NAMESPACE_STIR + +#endif diff --git a/src/include/stir/GammaImageFilter2D.h b/src/include/stir/GammaImageFilter2D.h new file mode 100644 index 0000000000..8113525e5d --- /dev/null +++ b/src/include/stir/GammaImageFilter2D.h @@ -0,0 +1,68 @@ +// +// +/*! + \file + \ingroup ImageProcessor + \brief Declaration of class stir::GammaImageFilter2D + + \author Dimitra Kyriakopoulou + \author Kris Thielemans +*/ +/* + Copyright (C) 2024, University College London + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +#ifndef __stir_GammaImageFilter2D_H__ +#define __stir_GammaImageFilter2D_H__ + +#include "stir/DataProcessor.h" +#include "stir/GammaArrayFilter2D.h" +#include "stir/DiscretisedDensity.h" +#include "stir/RegisteredParsingObject.h" + +START_NAMESPACE_STIR + +/*! + \ingroup ImageProcessor + \brief A class in the ImageProcessor hierarchy that implements gamma correction filtering. + + This class applies a gamma correction filter to a 3D volume using the GammaArrayFilter3D. + It is derived from RegisteredParsingObject, allowing it to be configured through parameter files. + + \tparam elemT The element type of the density (e.g., `float`, `double`). +*/ +template +class GammaImageFilter2D : public RegisteredParsingObject, + DataProcessor>, + DataProcessor>> +{ +private: + typedef RegisteredParsingObject, + DataProcessor>, + DataProcessor>> + base_type; + +public: + static const char* const registered_name; + + GammaImageFilter2D(); + +private: + GammaArrayFilter2D gamma_filter; + + void set_defaults() override; + void initialise_keymap() override; + + Succeeded virtual_set_up(const DiscretisedDensity<3, elemT>& density) override; + void virtual_apply(DiscretisedDensity<3, elemT>& density, const DiscretisedDensity<3, elemT>& in_density) const override; + void virtual_apply(DiscretisedDensity<3, elemT>& density) const override; +}; + +END_NAMESPACE_STIR + +#endif diff --git a/src/include/stir/ProjDataInMemory.h b/src/include/stir/ProjDataInMemory.h index 0ded4e017b..b101819b2f 100644 --- a/src/include/stir/ProjDataInMemory.h +++ b/src/include/stir/ProjDataInMemory.h @@ -62,6 +62,9 @@ class ProjDataInMemory : public ProjData //! Copy constructor ProjDataInMemory(const ProjDataInMemory& proj_data); + //! A static member to get the projection data in memory from a file + static shared_ptr read_from_file(const std::string& filename); + Viewgram get_viewgram(const int view_num, const int segment_num, const bool make_num_tangential_poss_odd = false, diff --git a/src/include/stir/Scanner.h b/src/include/stir/Scanner.h index 3fc5a6f26c..5b2210c694 100644 --- a/src/include/stir/Scanner.h +++ b/src/include/stir/Scanner.h @@ -165,6 +165,7 @@ class Scanner DiscoveryMI3ring, DiscoveryMI4ring, DiscoveryMI5ring, + DiscoveryMI6ring, HZLR, RATPET, PANDA, @@ -426,6 +427,7 @@ class Scanner inline float get_transaxial_block_spacing() const; /*! get total axial length covered by the detectors (incl. any gaps between blocks etc.) \todo Need to update this function when enabling different spacing between blocks and buckets etc. + For cylindrical scanners, calculates the length by the ring spacing. */ inline float get_axial_length() const; //@} (end of get block geometry info) diff --git a/src/include/stir/Scanner.inl b/src/include/stir/Scanner.inl index 2c742d21cc..050e80d172 100644 --- a/src/include/stir/Scanner.inl +++ b/src/include/stir/Scanner.inl @@ -298,6 +298,14 @@ Scanner::get_axial_block_spacing() const float Scanner::get_axial_length() const { + if (get_scanner_geometry() == "Cylindrical") + { + return get_num_rings() * get_ring_spacing(); + } + else if (get_scanner_geometry() == "Generic") + { + error("get_axial_length not yet implemented for Generic scanner geometry"); + } return get_num_axial_buckets() * get_num_axial_blocks_per_bucket() * get_axial_block_spacing() - (get_axial_block_spacing() - (get_num_axial_crystals_per_block() - 1) * get_axial_crystal_spacing()); } diff --git a/src/include/stir/WienerArrayFilter2D.h b/src/include/stir/WienerArrayFilter2D.h new file mode 100644 index 0000000000..624ecbbb0f --- /dev/null +++ b/src/include/stir/WienerArrayFilter2D.h @@ -0,0 +1,79 @@ +// +// +/*! + + \file + \ingroup Array + \brief Declaration of class stir::WienerArrayFilter2D + + \author Dimitra Kyriakopoulou + \author Kris Thielemans + +*/ +/* + Copyright (C) 2024, University College London + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +#ifndef __stir_WienerArrayFilter2D_H__ +#define __stir_WienerArrayFilter2D_H__ + +#include "stir/ArrayFunctionObject_2ArgumentImplementation.h" + +START_NAMESPACE_STIR + +/*! + \ingroup Array + \brief Applies a 2D Wiener filter on a 3D input array, slice by slice. + + This function applies a Wiener filter on each 2D slice of a 3D volume independently, using a fixed 3x3 window. + For each pixel in the slice, the filter estimates the local mean and variance, and uses these values, along with + the noise variance estimated from the entire slice, to reduce noise while preserving details. The filtered + output is stored in the `out_array`. + + The formula used for the Wiener filter is: + \f[ + \text{output}(i, j) = \left(\frac{\text{input}(i, j) - \text{localMean}(i, j)}{\max(\text{localVar}(i, j), \text{noise})}\right) + \cdot \max(\text{localVar}(i, j) - \text{noise}, 0) + \text{localMean}(i, j) + \f] + + - `localMean[i][j]` is the mean of the 3x3 neighborhood around pixel `(i, j)`. + - `localVar[i][j]` is the variance of the 3x3 neighborhood around pixel `(i, j)`. + - `noise` is the average noise variance estimated over the entire slice. + + + \warning The edges of each 2D slice are not processed, as the filter does not have sufficient neighboring pixels for the 3x3 + window. + + ### Wiener Filter Configuration + This filter is fully automated and does not require any parameters. + To enable it in the reconstruction process, include the following in the parameter file: + \code + post-filter type := Wiener + Wiener Filter Parameters := + End Wiener Filter Parameters := + \endcode + + \note This filter operates on each axial slice independently, and does not take into account neighboring slices. It is + effectively a 2D filter applied slice-by-slice on a 3D volume. + + */ + +template +class WienerArrayFilter2D : public ArrayFunctionObject_2ArgumentImplementation<3, elemT> +{ +public: + WienerArrayFilter2D(); + bool is_trivial() const override; + +private: + void do_it(Array<3, elemT>& out_array, const Array<3, elemT>& in_array) const override; +}; + +END_NAMESPACE_STIR + +#endif diff --git a/src/include/stir/WienerImageFilter2D.h b/src/include/stir/WienerImageFilter2D.h new file mode 100644 index 0000000000..36463362a3 --- /dev/null +++ b/src/include/stir/WienerImageFilter2D.h @@ -0,0 +1,69 @@ +// +// +/*! + + \file + \ingroup ImageProcessor + \brief Declaration of class stir::WienerImageFilter2D.h + + \author Dimitra Kyriakopoulou + \author Kris Thielemans + +*/ +/* + Copyright (C) 2024, University College London + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +#ifndef __stir_WienerImageFilter2D_H__ +# define __stir_WienerImageFilter2_H__ + +# include "stir/DataProcessor.h" +# include "stir/WienerArrayFilter2D.h" +# include "stir/DiscretisedDensity.h" +# include "stir/RegisteredParsingObject.h" + +START_NAMESPACE_STIR + +/*! + \ingroup ImageProcessor + \brief A class in the ImageProcessor hierarchy that implements wiener + filtering. + + As it is derived from RegisteredParsingObject, it implements all the + necessary things to parse parameter files etc. + */ +template +class WienerImageFilter2D : public RegisteredParsingObject, + DataProcessor>, + DataProcessor>> +{ +private: + typedef RegisteredParsingObject, + DataProcessor>, + DataProcessor>> + base_type; + +public: + static const char* const registered_name; + + WienerImageFilter2D(); + +private: + WienerArrayFilter2D wiener_filter; + + void set_defaults() override; + void initialise_keymap() override; + + Succeeded virtual_set_up(const DiscretisedDensity<3, elemT>& density) override; + void virtual_apply(DiscretisedDensity<3, elemT>& density, const DiscretisedDensity<3, elemT>& in_density) const override; + void virtual_apply(DiscretisedDensity<3, elemT>& density) const override; +}; + +END_NAMESPACE_STIR + +#endif diff --git a/src/include/stir/analytic/SRT2D/SRT2DReconstruction.h b/src/include/stir/analytic/SRT2D/SRT2DReconstruction.h new file mode 100644 index 0000000000..34d2699299 --- /dev/null +++ b/src/include/stir/analytic/SRT2D/SRT2DReconstruction.h @@ -0,0 +1,169 @@ +#ifndef __stir_analytic_SRT2D_SRT2DReconstruction_H__ +#define __stir_analytic_SRT2D_SRT2DReconstruction_H__ +/* + Copyright (C) 2024, University College London + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ +/*! + \file + \ingroup analytic + + \brief declares the stir::SRT2DReconstruction class + + \author Dimitra Kyriakopoulou + \author Kris Thielemans +*/ + +#include "stir/recon_buildblock/AnalyticReconstruction.h" +#include "stir/RegisteredParsingObject.h" +#include +#include +#include "stir/shared_ptr.h" + +START_NAMESPACE_STIR + +template +class DiscretisedDensity; +class Succeeded; +class ProjData; + +/*! \ingroup SRT2D + \brief Reconstruction class for 2D Spline Reconstruction Technique + + The reference for the implemented PET algorithm is: Fokas, A. S., A. Iserles, and V. Marinakis. "Reconstruction algorithm for +single photon emission computed tomography and its numerical implementation." *Journal of the Royal Society Interface* 3.6 (2006): +45-54. + + STIR implementations: Initial version June 2012, 1st updated version (4-point symmetry included) November 2012, 2nd updated +version (8-point symmetry included) July 2013, 3rd updated version 2014-2016, 4th updated version 2023-2024 + + \par Parameters + \verbatim +SRT2Dparameters := + +input file := input.hs +output filename prefix := output + +; output image parameters +; zoom defaults to 1 +zoom := -1 +; image size defaults to whole FOV +xy output image size (in pixels) := -1 + +; can be used to call SSRB first +; default means: call SSRB only if no axial compression is already present +;num segments to combine with ssrb := -1 + +END := + \endverbatim +*/ + +class SRT2DReconstruction + : public RegisteredParsingObject>, AnalyticReconstruction> +{ + // typedef AnalyticReconstruction base_type; + typedef RegisteredParsingObject>, AnalyticReconstruction> + base_type; +#ifdef SWIG + // work-around swig problem. It gets confused when using a private (or protected) + // typedef in a definition of a public typedef/member + public: +#else +private: +#endif + typedef DiscretisedDensity<3, float> TargetT; + +public: + //! Name which will be used when parsing a SRT2DReconstruction object + static const char* const registered_name; + + //! Default constructor (calls set_defaults()) + SRT2DReconstruction(); + /*! + \brief Constructor, initialises everything from parameter file, or (when + parameter_filename == "") by calling ask_parameters(). + */ + explicit SRT2DReconstruction(const std::string& parameter_filename); + + SRT2DReconstruction(const shared_ptr& proj_data_ptr_v, const int num_segments_to_combine = -1); + virtual std::string method_info() const; + + virtual void ask_parameters(); + + virtual Succeeded set_up(shared_ptr const& target_data_sptr); + +protected: // make parameters protected such that doc shows always up in doxygen + // parameters used for parsing + + //! number of segments to combine (with SSRB) before starting 2D reconstruction + /*! if -1, a value is chosen depending on the axial compression. + If there is no axial compression, num_segments_to_combine is + effectively set to 3, otherwise it is set to 1. + \see SSRB + */ + int num_segments_to_combine; + +private: + Succeeded actual_reconstruct(shared_ptr> const& target_image_ptr); + + virtual void set_defaults(); + virtual void initialise_keymap(); + virtual bool post_processing(); + + /*! + \brief Computes second derivatives for natural cubic spline interpolation. + + This function precomputes the second derivatives of the input data points + for use in cubic spline interpolation. The results are stored in the \a y2 vector. + + \param x Vector of x-coordinates of the input data points. + \param y Vector of y-coordinates of the input data points. + \param n The number of data points. + \param y2 Vector to store the computed second derivatives for spline interpolation. +*/ + void spline(const std::vector& x, const std::vector& y, int n, std::vector& y2) const; + + /*! + \brief Computes the Hilbert transform derivative for a set of projections. + + This function calculates the derivative of the Hilbert transform for a set of sampled data points. + It uses second derivatives, logarithmic differences, and a correction term to adjust the + computed derivative. + + \param x The x-coordinate for which the derivative is evaluated. + \param f Vector of function values at sampling points. + \param ddf Vector of second derivatives of the function. + \param p Vector of sampling positions. + \param sp The total number of sampling points. + \param lg Vector of logarithmic differences used for interpolation. + \param termC Correction term used for adjusting the derivative. + \return The computed Hilbert transform derivative at \a x. +*/ + float hilbert_der(float x, + const std::vector& f, + const std::vector& ddf, + const std::vector& p, + int sp, + const std::vector& lg, + float termC) const; + + /*! + \brief Performs numerical integration over a set of sampled data. + + This function uses a simple summation approach to numerically calculate integrals over tangential positions. + + \param dist The interval over which the integration is performed. + \param max The number of data points to integrate. + \param ff Vector containing the function values to be integrated. + \return The computed integral value. +*/ + float integ(float dist, int max, const std::vector& ff) const; +}; + +END_NAMESPACE_STIR + +#endif diff --git a/src/include/stir/analytic/SRT2DSPECT/SRT2DSPECTReconstruction.h b/src/include/stir/analytic/SRT2DSPECT/SRT2DSPECTReconstruction.h new file mode 100644 index 0000000000..240699fbc9 --- /dev/null +++ b/src/include/stir/analytic/SRT2DSPECT/SRT2DSPECTReconstruction.h @@ -0,0 +1,265 @@ +// +// +#ifndef __stir_analytic_SRT2DSPECT_SRT2DSPECTReconstruction_H__ +#define __stir_analytic_SRT2DSPECT_SRT2DSPECTReconstruction_H__ +/* + Copyright (C) 2024, University College London + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ +/*! + \file + \ingroup analytic + + \brief declares the stir::SRT2DSPECTReconstruction class + + \author Dimitra Kyriakopoulou + \author Kris Thielemans +*/ + +#include "stir/analytic/SRT2DSPECT/SRT2DSPECTReconstruction.h" +#include "stir/VoxelsOnCartesianGrid.h" +#include "stir/ProjDataInfoCylindricalArcCorr.h" +#include "stir/ArcCorrection.h" +#include "stir/SSRB.h" +#include "stir/ProjDataInMemory.h" +#include "stir/Array.h" +#include +#include "stir/Sinogram.h" +#include "stir/Viewgram.h" +#include +#include "stir/Bin.h" +#include "stir/round.h" +#include "stir/display.h" +#include +#include "stir/IO/interfile.h" +#include "stir/info.h" +#include + +#include "stir/recon_buildblock/AnalyticReconstruction.h" +#include "stir/RegisteredParsingObject.h" +#include +#include "stir/shared_ptr.h" + +START_NAMESPACE_STIR + +template +class DiscretisedDensity; +class Succeeded; +class ProjData; + +/*! \ingroup SRT2DSPECT + \brief Reconstruction class for 2D Spline Reconstruction Technique + + The reference for the implemented SPECT algorithm is: Fokas, A. S., A. Iserles, and V. Marinakis. "Reconstruction algorithm for +single photon emission computed tomography and its numerical implementation." *Journal of the Royal Society Interface* 3.6 (2006): +45-54. + + STIR implementations: initial version 2014-2016, 1st updated version 2023-2024 + + \par Parameters + + SRT2DSPECT takes two inputs: + - The emission sinogram, which represents the measured attenuated data. + - The attenuation projection sinogram, which is the Radon transform (line integrals) of the attenuation map. + + \verbatim +SRT2DSPECTparameters := + +input file := input.hs +attenuation projection filename := attenuation_projection_sinogram.hs +output filename prefix := output + +; output image parameters +; zoom defaults to 1 +zoom := -1 +; image size defaults to whole FOV +xy output image size (in pixels) := -1 + +; can be used to call SSRB first +; default means: call SSRB only if no axial compression is already present +;num segments to combine with ssrb := -1 + +END := + \endverbatim +*/ + +class SRT2DSPECTReconstruction : public RegisteredParsingObject>, + AnalyticReconstruction> +{ + // typedef AnalyticReconstruction base_type; + typedef RegisteredParsingObject>, AnalyticReconstruction> + base_type; +#ifdef SWIG + // work-around swig problem. It gets confused when using a private (or protected) + // typedef in a definition of a public typedef/member + public: +#else +private: +#endif + typedef DiscretisedDensity<3, float> TargetT; + +public: + //! Name which will be used when parsing a ProjectorByBinPair object + static const char* const registered_name; + + //! Default constructor (calls set_defaults()) + SRT2DSPECTReconstruction(); + /*! + \brief Constructor, initialises everything from parameter file, or (when + parameter_filename == "") by calling ask_parameters(). + */ + explicit SRT2DSPECTReconstruction(const std::string& parameter_filename); + + SRT2DSPECTReconstruction(const shared_ptr& proj_data_ptr_v, const int num_segments_to_combine = -1); + + virtual std::string method_info() const; + + virtual void ask_parameters(); + + virtual Succeeded set_up(shared_ptr const& target_data_sptr); + +protected: // make parameters protected such that doc shows always up in doxygen + // parameters used for parsing + + //! number of segments to combine (with SSRB) before starting 2D reconstruction + /*! if -1, a value is chosen depending on the axial compression. + If there is no axial compression, num_segments_to_combine is + effectively set to 3, otherwise it is set to 1. + \see SSRB + */ + int num_segments_to_combine; + std::string attenuation_projection_filename; + float thres_restr_bound; + std::vector thres_restr_bound_vector; + shared_ptr atten_data_ptr; + +private: + Succeeded actual_reconstruct(shared_ptr> const& target_image_ptr); + + virtual void set_defaults(); + virtual void initialise_keymap(); + virtual bool post_processing(); + + /*! + \brief Computes the Hilbert transform at a given node. + + This function calculates the Hilbert transform at a specific position + based on the given function values and their second derivatives. + + \param x The x-coordinate for the node. + \param f Vector of function values at sampling points. + \param ddf Vector of second derivatives of the function. + \param p Vector of sampling positions. + \param sp Number of sampling positions. + \param fn The function value at point \a x. + \return The computed Hilbert transform value at the specified node. + */ + float hilbert_node( + float x, const std::vector& f, const std::vector& ddf, const std::vector& p, int sp, float fn) const; + + /*! + \brief Computes the Hilbert transform for a set of sampled data. + + This function calculates the Hilbert transform for a single set of sampled data points, + which is used to adjust for phase shifts in projections during SPECT reconstruction. + The transform is computed using the function values, their second derivatives, + and logarithmic differences. + + \param x The x-coordinate where the transform is evaluated. + \param f Vector containing function values at the sampling positions. + \param ddf Vector of second derivatives of the function values. + \param p Vector of sampling positions. + \param sp The total number of sampling points. + \param lg Vector of logarithmic differences used for interpolation adjustments. + \return The computed Hilbert transform value at the specified x-coordinate. +*/ + float hilbert(float x, + const std::vector& f, + const std::vector& ddf, + const std::vector& p, + int sp, + std::vector& lg) const; + + /*! + \brief Computes the Hilbert transform derivatives for two sets of sampled data. + + This function calculates the Hilbert transform derivatives for two separate sets of + sampled data points, which are used to adjust for phase shifts in projections + during SPECT reconstruction. The derivatives are computed using the function values, + their second derivatives, and logarithmic differences. + + \param x The x-coordinate where the derivatives are evaluated. + \param f Vector containing function values for the first set. + \param ddf Vector of second derivatives for the first set. + \param f1 Vector containing function values for the second set. + \param ddf1 Vector of second derivatives for the second set. + \param p Vector of sampling positions. + \param sp Total number of sampling positions. + \param dhp Pointer to store the computed derivative for the first set. + \param dh1p Pointer to store the computed derivative for the second set. + \param lg Vector of logarithmic differences used for interpolation. + */ + void hilbert_der_double(float x, + const std::vector& f, + const std::vector& ddf, + const std::vector& f1, + const std::vector& ddf1, + const std::vector& p, + int sp, + float* dhp, + float* dh1p, + const std::vector& lg) const; + + /*! + \brief Performs interpolation using precomputed cubic splines. + + This function uses the precomputed second derivatives (calculated by the \a spline function) + to interpolate a value at a specified x-coordinate. It returns the interpolated value. + + \param xa The x-coordinates of the original data points. + \param ya The y-coordinates of the original data points. + \param y2a Precomputed second derivatives from the \a spline function. + \param n The number of data points. + \param x The x-coordinate where the interpolation is evaluated. + \return The interpolated function value at \a x. + + \note This function relies on the results of the \a spline function to efficiently compute the interpolation. +*/ + float splint(const std::vector& xa, const std::vector& ya, const std::vector& y2a, int n, float x) const; + + /*! + \brief Computes second derivatives for natural cubic spline interpolation. + + This function precomputes the second derivatives of the input data points + for use in cubic spline interpolation. The results are stored in the \a y2 vector. + + \param x The x-coordinates of the input data points. + \param y The y-coordinates of the input data points. + \param n The number of data points. + \param y2 Vector to store the computed second derivatives for spline interpolation. + + \note This function prepares the data for efficient use in the \a splint function. +*/ + void spline(const std::vector& x, const std::vector& y, int n, std::vector& y2) const; + + /*! + \brief Performs numerical integration over a specified range. + + This function computes numerical integration using a summation approach. + + \param dist The distance interval over which to perform the integration. + \param max The number of discrete data points used for the integration. + \param ff Array containing function values at each sampled point. + \return The computed integral over the specified range. + */ + float integ(float dist, int max, float ff[]) const; +}; + +END_NAMESPACE_STIR + +#endif diff --git a/src/include/stir/common.h b/src/include/stir/common.h index 5a501a040a..c28a0456d0 100644 --- a/src/include/stir/common.h +++ b/src/include/stir/common.h @@ -72,7 +72,9 @@ #include #include #include -#include +#if __cplusplus >= 202002L +# include +#endif //*************** namespace macros #define START_NAMESPACE_STIR \ @@ -128,10 +130,16 @@ //*************** START_NAMESPACE_STIR +#ifndef _PI +# if __cplusplus >= 202002L //! The constant pi to high precision. /*! \ingroup buildblock */ -#ifndef _PI -# define _PI boost::math::constants::pi() +# define _PI std::numbers::pi +# else +//! The constant pi to high precision. +/*! \ingroup buildblock */ +# define _PI 3.1415926535897932384626433832795028841971693993751058209749445923078164062 +# endif #endif //! Define the speed of light in mm / ps diff --git a/src/include/stir/interpolate_projdata.h b/src/include/stir/interpolate_projdata.h index 43889ca8a8..0fe97060c7 100644 --- a/src/include/stir/interpolate_projdata.h +++ b/src/include/stir/interpolate_projdata.h @@ -59,6 +59,8 @@ Succeeded interpolate_projdata(ProjData& proj_data_out, const ProjData& proj_data_in, const BasicCoordinate<3, BSpline::BSplineType>& these_types, const bool remove_interleaving); +Succeeded +interpolate_blocks_on_cylindrical_projdata(ProjData& proj_data_out, const ProjData& proj_data_in, bool remove_interleaving); //@} END_NAMESPACE_STIR diff --git a/src/include/stir/listmode/LmToProjData.h b/src/include/stir/listmode/LmToProjData.h index db858d7f31..5789d4e79d 100644 --- a/src/include/stir/listmode/LmToProjData.h +++ b/src/include/stir/listmode/LmToProjData.h @@ -202,6 +202,9 @@ class LmToProjData : public LmToProjDataAbstract bool get_store_prompts() const; void set_store_delayeds(bool); bool get_store_delayeds() const; + //! Returns the last processed timestamp in the listmode file + /*! This can be used to find the duration (in seconds) of the last time frame processed. */ + double get_last_processed_lm_rel_time() const; void set_num_segments_in_memory(int); int get_num_segments_in_memory() const; void set_num_events_to_store(long int); diff --git a/src/include/stir/modelling/PlasmaData.inl b/src/include/stir/modelling/PlasmaData.inl index e72f873309..ea74ff1995 100644 --- a/src/include/stir/modelling/PlasmaData.inl +++ b/src/include/stir/modelling/PlasmaData.inl @@ -21,6 +21,7 @@ #include "stir/warning.h" #include "stir/error.h" #include +#include START_NAMESPACE_STIR diff --git a/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl b/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl index 16f92a2e70..661744d503 100644 --- a/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl +++ b/src/include/stir/recon_buildblock/DataSymmetriesForBins_PET_CartesianGrid.inl @@ -216,6 +216,7 @@ DataSymmetriesForBins_PET_CartesianGrid::find_sym_op_bin0(int segment_num, int v // No symmetry is implemented for generic scanner return new TrivialSymmetryOperation(); } + return new TrivialSymmetryOperation(); } // from symmetries @@ -390,6 +391,8 @@ DataSymmetriesForBins_PET_CartesianGrid::find_sym_op_general_bin(int s, int segm // No symmetry is implemented for generic scanner return new TrivialSymmetryOperation(); } + + return new TrivialSymmetryOperation(); } bool diff --git a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h index f7fa35ce6e..9fcbc0650e 100644 --- a/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h +++ b/src/include/stir/recon_buildblock/ML_estimate_component_based_normalisation.h @@ -9,10 +9,11 @@ /*! \file \ingroup recon_buildblock - \brief Declaration of ML_estimate_component_based_normalisation + \brief Declaration of stir::ML_estimate_component_based_normalisation \author Kris Thielemans */ #include "stir/common.h" +#include START_NAMESPACE_STIR diff --git a/src/include/stir/scatter/ScatterEstimation.h b/src/include/stir/scatter/ScatterEstimation.h index 38df94df48..84d91e745c 100644 --- a/src/include/stir/scatter/ScatterEstimation.h +++ b/src/include/stir/scatter/ScatterEstimation.h @@ -39,6 +39,7 @@ #include "stir/stir_math.h" #include "stir/FilePath.h" +#include "stir/recon_buildblock/ForwardProjectorByBin.h" START_NAMESPACE_STIR @@ -119,7 +120,8 @@ class ScatterEstimation : public ParsingObject //! make projdata 2D shared pointer shared_ptr make_2D_projdata_sptr(const shared_ptr in_3d_sptr); - shared_ptr make_2D_projdata_sptr(const shared_ptr in_3d_sptr, string template_filename); + shared_ptr + make_2D_projdata_sptr(const shared_ptr in_3d_sptr, string template_filename, const bool do_normalisation = false); //! //! \brief set_up @@ -170,10 +172,21 @@ class ScatterEstimation : public ParsingObject void set_recompute_mask_image(bool arg); void set_recompute_mask_projdata(bool arg); + //! Sets the forward projector for the scatter estimation + void set_forward_projector_for_mask_sptr(const shared_ptr projector_sptr); + //! Get the forward projector used for the scatter estimation mask calculation + shared_ptr get_forward_projector_for_mask_sptr() const; + inline void set_scatter_simulation_method_sptr(const shared_ptr); inline void set_num_iterations(int); + inline unsigned int get_half_filter_width() const; + inline void set_half_filter_width(unsigned int); + + inline void + set_downsample_scanner(bool downsample_scanner, int downsampled_number_of_rings = -1, int downsampled_detectors_per_ring = -1); + void set_output_scatter_estimate_prefix(const std::string&); void set_export_scatter_estimates_of_each_iteration(bool); @@ -269,6 +282,8 @@ class ScatterEstimation : public ParsingObject //! The file name for the attenuation coefficients. //! If they are to be recalculated they will be stored here, if set. std::string atten_coeff_filename; + //! ForwardProjector + shared_ptr forward_projector_for_mask_sptr; //! \details the set of parameters to obtain a mask from the attenuation image /*! The attenuation image will be thresholded to find a plausible mask for where there @@ -382,6 +397,9 @@ class ScatterEstimation : public ParsingObject float min_scale_value; bool downsample_scanner_bool; + int downsampled_number_of_rings; + int downsampled_detectors_per_ring; + //! unsigned int half_filter_width; diff --git a/src/include/stir/scatter/ScatterEstimation.inl b/src/include/stir/scatter/ScatterEstimation.inl index 73336a24d8..dc8bd9f2bd 100644 --- a/src/include/stir/scatter/ScatterEstimation.inl +++ b/src/include/stir/scatter/ScatterEstimation.inl @@ -66,4 +66,26 @@ ScatterEstimation::set_num_iterations(int arg) this->num_scatter_iterations = arg; } +unsigned int +ScatterEstimation::get_half_filter_width() const +{ + return this->half_filter_width; +} + +void +ScatterEstimation::set_half_filter_width(unsigned int arg) +{ + this->half_filter_width = arg; +} + +void +ScatterEstimation::set_downsample_scanner(bool downsample_scanner, + int downsampled_number_of_rings, + int downsampled_detectors_per_ring) +{ + this->downsample_scanner_bool = downsample_scanner; + this->downsampled_number_of_rings = downsampled_number_of_rings; + this->downsampled_detectors_per_ring = downsampled_detectors_per_ring; +} + END_NAMESPACE_STIR diff --git a/src/include/stir/warning.h b/src/include/stir/warning.h index 7d7fc56787..317aee04b7 100644 --- a/src/include/stir/warning.h +++ b/src/include/stir/warning.h @@ -17,11 +17,9 @@ \author Kris Thielemans */ -#include "stir/common.h" #include "stir/Verbosity.h" -#include - -#include "TextWriter.h" +#include "stir/TextWriter.h" +#include START_NAMESPACE_STIR diff --git a/src/listmode_buildblock/LmToProjData.cxx b/src/listmode_buildblock/LmToProjData.cxx index edcce9519b..25dc99b898 100644 --- a/src/listmode_buildblock/LmToProjData.cxx +++ b/src/listmode_buildblock/LmToProjData.cxx @@ -184,6 +184,12 @@ LmToProjData::get_output_filename_prefix() const return output_filename_prefix; } +double +LmToProjData::get_last_processed_lm_rel_time() const +{ + return this->current_time; +} + void LmToProjData::set_output_projdata_sptr(shared_ptr& arg) { @@ -442,8 +448,15 @@ LmToProjData::set_up() else if (frame_defs.get_num_frames() < 1) { // make a single frame starting from 0. End value will be ignored. - vector> frame_times(1, pair(0, 0)); - frame_defs = TimeFrameDefinitions(frame_times); + frame_defs = lm_data_ptr->get_exam_info_sptr()->get_time_frame_definitions(); + if (frame_defs.get_num_time_frames() == 0) + { + vector> frame_times(1, pair(0, 0)); + frame_defs = TimeFrameDefinitions(frame_times); + } + if (num_events_to_store != 0) + warning( + "LmToProjData: num_events_to_store has been selected. The frame duration in the Interfile header will be incorrect!"); } return Succeeded::yes; diff --git a/src/recon_buildblock/BinNormalisationSPECT.cxx b/src/recon_buildblock/BinNormalisationSPECT.cxx index 43971d7159..5b3ec1689e 100644 --- a/src/recon_buildblock/BinNormalisationSPECT.cxx +++ b/src/recon_buildblock/BinNormalisationSPECT.cxx @@ -318,7 +318,7 @@ BinNormalisationSPECT::read_uniformity_table(Array<3, float>& uniformity) const for (int n = 1; n <= num_detector_heads; n++) { - const std::string n_string = boost::lexical_cast(n); + const std::string n_string = std::to_string(n); const std::string filename(this->folder_prefix + n_string + "/" + uniformity_filename); std::ifstream input(filename.c_str()); diff --git a/src/recon_buildblock/CMakeLists.txt b/src/recon_buildblock/CMakeLists.txt index 61286275f7..ae64f14fce 100644 --- a/src/recon_buildblock/CMakeLists.txt +++ b/src/recon_buildblock/CMakeLists.txt @@ -147,19 +147,19 @@ endif() target_link_libraries(recon_buildblock PUBLIC modelling_buildblock display numerics_buildblock listmode_buildblock data_buildblock buildblock spatial_transformation_buildblock ) if (STIR_WITH_NiftyPET_PROJECTOR) - target_link_libraries(recon_buildblock PUBLIC NiftyPET::petprj) - target_link_libraries(recon_buildblock PUBLIC NiftyPET::mmr_auxe) - target_link_libraries(recon_buildblock PUBLIC NiftyPET::mmr_lmproc) - target_link_libraries(recon_buildblock PUBLIC CUDA::cudart) + target_link_libraries(recon_buildblock PRIVATE NiftyPET::petprj) + target_link_libraries(recon_buildblock PRIVATE NiftyPET::mmr_auxe) + target_link_libraries(recon_buildblock PRIVATE NiftyPET::mmr_lmproc) + target_link_libraries(recon_buildblock PRIVATE CUDA::cudart) endif() if (STIR_WITH_Parallelproj_PROJECTOR) - target_link_libraries(recon_buildblock PUBLIC parallelproj::parallelproj_c) + target_link_libraries(recon_buildblock PRIVATE parallelproj::parallelproj_c) if (parallelproj_built_with_CUDA) - target_link_libraries(recon_buildblock PUBLIC parallelproj::parallelproj_cuda) + target_link_libraries(recon_buildblock PRIVATE parallelproj::parallelproj_cuda) endif() endif() if (STIR_WITH_CUDA) - target_link_libraries(recon_buildblock PUBLIC CUDA::cudart) -endif() \ No newline at end of file + target_link_libraries(recon_buildblock PRIVATE CUDA::cudart) +endif() diff --git a/src/scatter_buildblock/CreateTailMaskFromACFs.cxx b/src/scatter_buildblock/CreateTailMaskFromACFs.cxx index 65a249dc0f..884b83a5c2 100644 --- a/src/scatter_buildblock/CreateTailMaskFromACFs.cxx +++ b/src/scatter_buildblock/CreateTailMaskFromACFs.cxx @@ -11,6 +11,7 @@ #include "stir/scatter/CreateTailMaskFromACFs.h" #include "stir/Bin.h" #include "boost/lambda/lambda.hpp" +#include "boost/format.hpp" #include "stir/error.h" START_NAMESPACE_STIR @@ -147,8 +148,8 @@ CreateTailMaskFromACFs::process_data() count += mask_left_size + mask_right_size; #endif } - std::cout << count << " bins in mask for sinogram at segment " << bin.segment_num() << ", axial_pos " - << bin.axial_pos_num() << "\n"; + info(boost::format("%1% bins in mask for sinogram at segment %2%, axial_pos %3%") % count % bin.segment_num() + % bin.axial_pos_num()); if (this->mask_proj_data->set_sinogram(mask_sinogram) != Succeeded::yes) return Succeeded::no; } diff --git a/src/scatter_buildblock/ScatterEstimation.cxx b/src/scatter_buildblock/ScatterEstimation.cxx index c157f4af81..54cf56167e 100644 --- a/src/scatter_buildblock/ScatterEstimation.cxx +++ b/src/scatter_buildblock/ScatterEstimation.cxx @@ -23,6 +23,7 @@ #include "stir/recon_buildblock/ChainedBinNormalisation.h" #include "stir/ProjDataInterfile.h" #include "stir/ProjDataInMemory.h" +#include "stir/inverse_SSRB.h" #include "stir/ExamInfo.h" #include "stir/ProjDataInfo.h" #include "stir/ProjDataInfoCylindricalNoArcCorr.h" @@ -77,6 +78,8 @@ ScatterEstimation::set_defaults() this->override_scanner_template = true; this->override_density_image = true; this->downsample_scanner_bool = true; + this->downsampled_number_of_rings = -1; + this->downsampled_detectors_per_ring = -1; this->remove_interleaving = true; this->atten_image_filename = ""; this->atten_coeff_filename = ""; @@ -109,6 +112,10 @@ ScatterEstimation::initialise_keymap() this->parser.add_key("mask projdata filename", &this->mask_projdata_filename); this->parser.add_key("tail fitting parameter filename", &this->tail_mask_par_filename); // END MASK + + // Forward projector for mask projection + this->parser.add_parsing_key("forward projector for mask type", &this->forward_projector_for_mask_sptr); + this->parser.add_key("background projdata filename", &this->back_projdata_filename); this->parser.add_parsing_key("Normalisation type", &this->norm_3d_sptr); this->parser.add_key("attenuation correction factors filename", &this->atten_coeff_filename); @@ -124,6 +131,8 @@ ScatterEstimation::initialise_keymap() this->parser.add_parsing_key("Scatter Simulation type", &this->scatter_simulation_sptr); this->parser.add_key("scatter simulation parameter filename", &this->scatter_sim_par_filename); this->parser.add_key("use scanner downsampling in scatter simulation", &this->downsample_scanner_bool); + this->parser.add_key("override number of downsampled rings", &this->downsampled_number_of_rings); + this->parser.add_key("override number of downsampled detectors per ring", &this->downsampled_detectors_per_ring); this->parser.add_key("override attenuation image", &this->override_density_image); this->parser.add_key("override scanner template", &this->override_scanner_template); @@ -182,7 +191,9 @@ ScatterEstimation::make_2D_projdata_sptr(const shared_ptr in_3d_sptr) } shared_ptr -ScatterEstimation::make_2D_projdata_sptr(const shared_ptr in_3d_sptr, string template_filename) +ScatterEstimation::make_2D_projdata_sptr(const shared_ptr in_3d_sptr, + string template_filename, + const bool do_normalisation) { shared_ptr out_2d_sptr; if (in_3d_sptr->get_proj_data_info_sptr()->get_scanner_sptr()->get_scanner_geometry() == "Cylindrical") @@ -193,7 +204,7 @@ ScatterEstimation::make_2D_projdata_sptr(const shared_ptr in_3d_sptr, this->input_projdata_2d_sptr->get_exam_info_sptr(), this->input_projdata_2d_sptr->get_proj_data_info_sptr()->create_shared_clone()); - SSRB(*out_2d_sptr, *in_3d_sptr, false); + SSRB(*out_2d_sptr, *in_3d_sptr, do_normalisation); } else { @@ -444,6 +455,18 @@ ScatterEstimation::set_recompute_mask_projdata(bool arg) this->recompute_mask_projdata = arg; } +void +ScatterEstimation::set_forward_projector_for_mask_sptr(const shared_ptr projector_sptr) +{ + this->forward_projector_for_mask_sptr = projector_sptr; +} + +shared_ptr +ScatterEstimation::get_forward_projector_for_mask_sptr() const +{ + return forward_projector_for_mask_sptr; +} + bool ScatterEstimation::already_setup() const { @@ -597,7 +620,7 @@ ScatterEstimation::set_up() } if (this->downsample_scanner_bool) - this->scatter_simulation_sptr->downsample_scanner(); + this->scatter_simulation_sptr->downsample_scanner(this->downsampled_number_of_rings, this->downsampled_detectors_per_ring); // Check if Load a mask proj_data @@ -653,7 +676,7 @@ ScatterEstimation::set_up_iterative(shared_ptrmultiplicative_binnorm_2d_sptr.reset(new ChainedBinNormalisation(norm_coeff_2d_sptr, atten_coeff_2d_sptr)); this->multiplicative_binnorm_2d_sptr->set_up( - this->back_projdata_sptr->get_exam_info_sptr(), + this->input_projdata_2d_sptr->get_exam_info_sptr(), this->input_projdata_2d_sptr->get_proj_data_info_sptr()->create_shared_clone()); iterative_object->get_objective_function_sptr()->set_normalisation_sptr(multiplicative_binnorm_2d_sptr); } @@ -834,7 +857,7 @@ ScatterEstimation::process_data() float local_min_scale_value = 0.5f; float local_max_scale_value = 0.5f; - stir::BSpline::BSplineType spline_type = stir::BSpline::quadratic; + stir::BSpline::BSplineType spline_type = stir::BSpline::linear; // This has been set to 2D or 3D in the set_up() shared_ptr unscaled_est_projdata_sptr( @@ -1026,16 +1049,13 @@ ScatterEstimation::process_data() shared_ptr normalisation_factors_3d_sptr = this->get_normalisation_object_sptr(this->multiplicative_binnorm_sptr); - upsample_and_fit_scatter_estimate(*scatter_estimate_sptr, - *this->input_projdata_sptr, - *temp_projdata, - *normalisation_factors_3d_sptr, - *this->input_projdata_sptr, - 1.0f, - 1.0f, - 1, - spline_type, - false); + ProjDataInMemory interpolated_scatter(this->input_projdata_sptr->get_exam_info_sptr(), + this->input_projdata_sptr->get_proj_data_info_sptr()->create_shared_clone()); + inverse_SSRB(interpolated_scatter, *temp_projdata); + normalisation_factors_3d_sptr->set_up(this->input_projdata_sptr->get_exam_info_sptr(), + this->input_projdata_sptr->get_proj_data_info_sptr()->create_shared_clone()); + normalisation_factors_3d_sptr->undo(interpolated_scatter); + scatter_estimate_sptr->fill(interpolated_scatter); } else { @@ -1233,30 +1253,32 @@ ScatterEstimation::project_mask_image() } } - shared_ptr forw_projector_sptr; - shared_ptr PM(new ProjMatrixByBinUsingRayTracing()); - forw_projector_sptr.reset(new ForwardProjectorByBinUsingProjMatrixByBin(PM)); + if (!this->forward_projector_for_mask_sptr) + { + shared_ptr PM(new ProjMatrixByBinUsingRayTracing()); + forward_projector_for_mask_sptr.reset(new ForwardProjectorByBinUsingProjMatrixByBin(PM)); + } info(boost::format("ScatterEstimation: Forward projector used for the calculation of " "the tail mask: %1%") - % forw_projector_sptr->parameter_info()); + % forward_projector_for_mask_sptr->parameter_info()); shared_ptr mask_projdata; if (run_in_2d_projdata) { - forw_projector_sptr->set_up(this->input_projdata_2d_sptr->get_proj_data_info_sptr(), this->mask_image_sptr); + forward_projector_for_mask_sptr->set_up(this->input_projdata_2d_sptr->get_proj_data_info_sptr(), this->mask_image_sptr); mask_projdata.reset(new ProjDataInMemory(this->input_projdata_2d_sptr->get_exam_info_sptr(), this->input_projdata_2d_sptr->get_proj_data_info_sptr())); } else { - forw_projector_sptr->set_up(this->input_projdata_sptr->get_proj_data_info_sptr(), this->mask_image_sptr); + forward_projector_for_mask_sptr->set_up(this->input_projdata_sptr->get_proj_data_info_sptr(), this->mask_image_sptr); mask_projdata.reset(new ProjDataInMemory(this->input_projdata_sptr->get_exam_info_sptr(), this->input_projdata_sptr->get_proj_data_info_sptr())); } - forw_projector_sptr->forward_project(*mask_projdata, *this->mask_image_sptr); + forward_projector_for_mask_sptr->forward_project(*mask_projdata, *this->mask_image_sptr); // add 1 to be able to use create_tail_mask_from_ACFs (which expects ACFs, // so complains if the threshold is too low) diff --git a/src/scatter_buildblock/ScatterSimulation.cxx b/src/scatter_buildblock/ScatterSimulation.cxx index 8411b3f53c..522178b83a 100644 --- a/src/scatter_buildblock/ScatterSimulation.cxx +++ b/src/scatter_buildblock/ScatterSimulation.cxx @@ -841,19 +841,51 @@ ScatterSimulation::downsample_scanner(int new_num_rings, int new_num_dets) float approx_num_non_arccorrected_bins; if (new_scanner_sptr->get_scanner_geometry() != "Cylindrical") { - new_num_dets = this->proj_data_info_sptr->get_scanner_ptr()->get_num_detectors_per_ring(); - approx_num_non_arccorrected_bins = this->proj_data_info_sptr->get_num_tangential_poss(); + if (new_num_dets <= 0) + { // by default, do not downsample the detectors per ring for BlocksOnCylindrical + new_num_dets = this->proj_data_info_sptr->get_scanner_ptr()->get_num_detectors_per_ring(); + } + // STIR does not like block spacings that are smaller than number of crystals times crystal spacing, + // therefore add a 1% extension on top for the downsampled scanner, to avoid running into floating point issues + const float block_spacing_factor = 1.01; + + // extend the bins by a small amount to avoid edge-effects, at the expense of longer computation time + approx_num_non_arccorrected_bins = ceil(this->proj_data_info_sptr->get_num_tangential_poss() * float(new_num_dets) + / old_scanner_ptr->get_num_detectors_per_ring()) + + 1; // preserve the length of the scanner, but place crystals equidistantly float scanner_length = old_scanner_ptr->get_axial_length(); float new_ring_spacing = scanner_length / (new_num_rings - 1); - new_scanner_sptr->set_num_axial_blocks_per_bucket(1); + new_scanner_sptr->set_num_transaxial_blocks_per_bucket(1); + new_scanner_sptr->set_num_rings(new_num_rings); + new_scanner_sptr->set_num_detectors_per_ring(new_num_dets); + new_scanner_sptr->set_axial_crystal_spacing(new_ring_spacing); new_scanner_sptr->set_ring_spacing(new_ring_spacing); new_scanner_sptr->set_num_axial_crystals_per_block(new_num_rings); - new_scanner_sptr->set_axial_block_spacing(new_ring_spacing * new_scanner_sptr->get_num_axial_crystals_per_block()); + new_scanner_sptr->set_axial_block_spacing(new_ring_spacing * new_num_rings * block_spacing_factor); + + float transaxial_bucket_width + = old_scanner_ptr->get_transaxial_block_spacing() * (old_scanner_ptr->get_num_transaxial_blocks_per_bucket() - 1) + + old_scanner_ptr->get_transaxial_crystal_spacing() * (old_scanner_ptr->get_num_transaxial_crystals_per_block() - 1); + + int num_trans_buckets = old_scanner_ptr->get_num_transaxial_buckets(); + // get a new number of detectors that is a multiple of the number of buckets to preserve scanner shape + new_scanner_sptr->set_num_detectors_per_ring(new_num_dets); + int new_transaxial_dets_per_bucket = new_num_dets / num_trans_buckets; + float new_det_spacing = transaxial_bucket_width / (new_transaxial_dets_per_bucket - 1); + + new_scanner_sptr->set_num_transaxial_blocks_per_bucket(1); + new_scanner_sptr->set_num_transaxial_crystals_per_block(new_transaxial_dets_per_bucket); + new_scanner_sptr->set_transaxial_crystal_spacing(new_det_spacing); + new_scanner_sptr->set_transaxial_block_spacing(new_transaxial_dets_per_bucket * new_det_spacing * block_spacing_factor); + // avoid problems with Scanner checks by setting singles_units to 1 crystal + // (only used for dead-time processing in ECAY norm) + new_scanner_sptr->set_num_axial_crystals_per_singles_unit(1); + new_scanner_sptr->set_num_transaxial_crystals_per_singles_unit(1); } else { diff --git a/src/test/CMakeLists.txt b/src/test/CMakeLists.txt index 6b7ae72fda..e111be65b0 100644 --- a/src/test/CMakeLists.txt +++ b/src/test/CMakeLists.txt @@ -74,7 +74,7 @@ set(buildblock_simple_tests include(stir_test_exe_targets) foreach(source ${buildblock_simple_tests}) - create_stir_test(${source} "buildblock;IO;buildblock;numerics_buildblock;display;IO;recon_buildblock;Shape_buildblock" "") + create_stir_test(${source} "buildblock;IO;buildblock;numerics_buildblock;display;IO;recon_buildblock;Shape_buildblock;scatter_buildblock" "") endforeach() # diff --git a/src/test/test_interpolate_projdata.cxx b/src/test/test_interpolate_projdata.cxx index 29c554ce7d..3153678b3e 100644 --- a/src/test/test_interpolate_projdata.cxx +++ b/src/test/test_interpolate_projdata.cxx @@ -41,6 +41,7 @@ #include "stir/Shape/Box3D.h" #include "stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h" #include "stir/recon_buildblock/ForwardProjectorByBinUsingProjMatrixByBin.h" +#include "stir/scatter/SingleScatterSimulation.h" #include "stir/RunTests.h" @@ -56,6 +57,8 @@ class InterpolationTests : public RunTests void scatter_interpolation_test_cyl(); void scatter_interpolation_test_blocks_asymmetric(); void scatter_interpolation_test_cyl_asymmetric(); + void scatter_interpolation_test_blocks_downsampled(); + void transaxial_upsampling_interpolation_test_blocks(); void check_symmetry(const SegmentBySinogram& segment); void compare_segment(const SegmentBySinogram& segment1, const SegmentBySinogram& segment2, float maxDiff); @@ -675,6 +678,280 @@ InterpolationTests::scatter_interpolation_test_cyl_asymmetric() compare_segment_shape(full_size_model_sino.get_segment_by_sinogram(0), interpolated_proj_data.get_segment_by_sinogram(0), 2); } +void +InterpolationTests::scatter_interpolation_test_blocks_downsampled() +{ + info("Performing downampled interpolation test for BlocksOnCylindrical scanner"); + auto time_frame_def = TimeFrameDefinitions(); + time_frame_def.set_num_time_frames(1); + time_frame_def.set_time_frame(1, 0, 1e9); + auto exam_info = ExamInfo(); + exam_info.set_high_energy_thres(650); + exam_info.set_low_energy_thres(425); + exam_info.set_time_frame_definitions(time_frame_def); + + // define the original scanner and a downsampled one, as it would be used for scatter simulation + auto scanner = Scanner(Scanner::User_defined_scanner, + "Some_BlocksOnCylindrical_Scanner", + 96, + 30, + int(150 * 96 / 192), + int(150 * 96 / 192), + 127, + 6.5, + 3.313, + 4.156, + -3.1091819, + 5, + 3, + 6, + 4, + 1, + 1, + 1, + 0.14, + 511, + 1, + 0, + 500, + "BlocksOnCylindrical", + 3.313, + 7.0, + 20.0, + 29.0); + auto downsampled_scanner = Scanner(Scanner::User_defined_scanner, + "Some_Downsampled_BlocksOnCylindrical_Scanner", + 64, + 8, + int(150 * 64 / 192 + 1), + int(150 * 64 / 192 + 1), + 127, + 6.5, + 16.652, + 6.234, + -3.1091819, + 1, + 1, + 8, + 8, + 1, + 1, + 1, + 0.17, + 511, + 1, + 0, + 500, + "BlocksOnCylindrical", + 13.795, + 15.4286, + 112.0, + 125.0); + + auto proj_data_info = shared_ptr(std::move( + ProjDataInfo::construct_proj_data_info(std::make_shared(scanner), 1, 29, 48, int(150 * 96 / 192), false))); + auto downsampled_proj_data_info = shared_ptr(std::move(ProjDataInfo::construct_proj_data_info( + std::make_shared(downsampled_scanner), 1, 0, 32, int(150 * 64 / 192 + 1), false))); + + auto proj_data = ProjDataInMemory(std::make_shared(exam_info), proj_data_info); + auto downsampled_proj_data = ProjDataInMemory(std::make_shared(exam_info), downsampled_proj_data_info); + + // define a cylinder and a box that are off-centre, such that the shapes in the sinogram can be compared + auto emission_map = VoxelsOnCartesianGrid(*proj_data_info, 1); + auto cyl_map = VoxelsOnCartesianGrid(*proj_data_info, 1); + auto cylinder = EllipsoidalCylinder(40, 40, 20, CartesianCoordinate3D(80, 100, 0)); + cylinder.construct_volume(cyl_map, CartesianCoordinate3D(1, 1, 1)); + auto box = Box3D(20, 20, 20, CartesianCoordinate3D(30, -20, 70)); + box.construct_volume(emission_map, CartesianCoordinate3D(1, 1, 1)); + emission_map += cyl_map; + + // project the emission map onto the full-scale scanner proj data + auto pm = ProjMatrixByBinUsingRayTracing(); + pm.set_use_actual_detector_boundaries(true); + pm.enable_cache(false); + auto forw_proj = ForwardProjectorByBinUsingProjMatrixByBin(std::make_shared(pm)); + forw_proj.set_up(proj_data_info, std::make_shared>(emission_map)); + auto full_size_model_sino = ProjDataInMemory(proj_data); + full_size_model_sino.fill(0); + forw_proj.forward_project(full_size_model_sino, emission_map); + + // also project onto the downsampled scanner + emission_map = VoxelsOnCartesianGrid(*downsampled_proj_data_info, 1); + cyl_map = VoxelsOnCartesianGrid(*downsampled_proj_data_info, 1); + cylinder.construct_volume(cyl_map, CartesianCoordinate3D(1, 1, 1)); + box.construct_volume(emission_map, CartesianCoordinate3D(1, 1, 1)); + emission_map += cyl_map; + forw_proj.set_up(downsampled_proj_data_info, std::make_shared>(emission_map)); + auto downsampled_model_sino = ProjDataInMemory(downsampled_proj_data); + downsampled_model_sino.fill(0); + forw_proj.forward_project(downsampled_model_sino, emission_map); + + // write the proj data to file + downsampled_model_sino.write_to_file("transaxially_downsampled_sino.hs"); + full_size_model_sino.write_to_file("transaxially_full_size_sino.hs"); + + // interpolate the downsampled proj data to the original scanner size and fill in oblique sinograms + auto interpolated_direct_proj_data = ProjDataInMemory(proj_data); + interpolate_projdata(interpolated_direct_proj_data, downsampled_model_sino, BSpline::linear, false); + auto interpolated_proj_data = ProjDataInMemory(proj_data); + inverse_SSRB(interpolated_proj_data, interpolated_direct_proj_data); + + // write the proj data to file + interpolated_proj_data.write_to_file("transaxially_interpolated_sino.hs"); + + // compare to ground truth + compare_segment_shape(full_size_model_sino.get_segment_by_sinogram(0), interpolated_proj_data.get_segment_by_sinogram(0), 3); +} + +void +InterpolationTests::transaxial_upsampling_interpolation_test_blocks() +{ + info("Performing transaxial downampled interpolation test for BlocksOnCylindrical scanner"); + auto time_frame_def = TimeFrameDefinitions(); + time_frame_def.set_num_time_frames(1); + time_frame_def.set_time_frame(1, 0, 1e9); + auto exam_info = ExamInfo(); + exam_info.set_high_energy_thres(650); + exam_info.set_low_energy_thres(425); + exam_info.set_time_frame_definitions(time_frame_def); + + // define the original scanner and a downsampled one, as it would be used for scatter simulation + auto scanner = Scanner(Scanner::User_defined_scanner, + "Some_BlocksOnCylindrical_Scanner", + 96, + 3, + 60, + 60, + 127, + 6.5, + 3.313, + 1.65, + -3.1091819, + 1, + 3, + 3, + 4, + 1, + 1, + 1, + 0.14, + 511, + 1, + 0, + 500, + "BlocksOnCylindrical", + 3.313, + 7.0, + 20.0, + 29.0); + auto proj_data_info = shared_ptr( + std::move(ProjDataInfo::construct_proj_data_info(std::make_shared(scanner), 1, 0, 48, 60, false))); + + // use the code in scatter simulation to downsample the scanner + auto scatter_simulation = SingleScatterSimulation(); + scatter_simulation.set_template_proj_data_info(*proj_data_info); + scatter_simulation.set_exam_info(exam_info); + scatter_simulation.downsample_scanner(-1, 96 / 4); // number of detectors per ring reduced by factor of four + auto downsampled_proj_data_info = scatter_simulation.get_template_proj_data_info_sptr(); + + auto proj_data = ProjDataInMemory(std::make_shared(exam_info), proj_data_info); + auto downsampled_proj_data = ProjDataInMemory(std::make_shared(exam_info), downsampled_proj_data_info); + + // define a cylinder precisely in the middle of the FOV + auto emission_map = VoxelsOnCartesianGrid(*downsampled_proj_data_info, 1); + make_symmetric_object(emission_map); + + // project the cylinder onto the full-scale scanner proj data + auto pm = ProjMatrixByBinUsingRayTracing(); + pm.set_use_actual_detector_boundaries(true); + pm.enable_cache(false); + auto forw_proj = ForwardProjectorByBinUsingProjMatrixByBin(std::make_shared(pm)); + forw_proj.set_up(proj_data_info, std::make_shared>(emission_map)); + auto full_size_model_sino = ProjDataInMemory(proj_data); + full_size_model_sino.fill(0); + forw_proj.forward_project(full_size_model_sino, emission_map); + + // also project onto the downsampled scanner + emission_map = VoxelsOnCartesianGrid(*downsampled_proj_data_info, 1); + make_symmetric_object(emission_map); + forw_proj.set_up(downsampled_proj_data_info, std::make_shared>(emission_map)); + auto downsampled_model_sino = ProjDataInMemory(downsampled_proj_data); + downsampled_model_sino.fill(0); + forw_proj.forward_project(downsampled_model_sino, emission_map); + + // write the proj data to file + downsampled_model_sino.write_to_file("transaxially_downsampled_sino_for_LOR.hs"); + + // interpolate the downsampled proj data to the original scanner size and fill in oblique sinograms + auto interpolated_direct_proj_data = ProjDataInMemory(proj_data); + interpolate_projdata(interpolated_direct_proj_data, downsampled_model_sino, BSpline::linear, false); + auto interpolated_proj_data = ProjDataInMemory(proj_data); + inverse_SSRB(interpolated_proj_data, interpolated_direct_proj_data); + + // write the proj data to file + interpolated_proj_data.write_to_file("transaxially_interpolated_sino_for_LOR.hs"); + + // Identify the bins which should be identical between the downsampled and the interpolated sinogram: + // Each module has 96 / 8 = 12 crystal in the full size scanner, organised in 3 blocks of 4 crystals, while + // the downsampled scanner has 3 crystals per module. The idea is that the centre of the outer two + // is in exactly the same position than the centre of the first and last crystal in the full size scanner. + SegmentBySinogram sinogram_downsampled = downsampled_proj_data.get_empty_segment_by_sinogram(0, false, 0); + SegmentBySinogram sinogram_full_size = proj_data.get_empty_segment_by_sinogram(0, false, 0); + const auto pdi_downsampled + = dynamic_cast(&(*downsampled_proj_data.get_proj_data_info_sptr())); + const auto pdi_full_size = dynamic_cast(&(*sinogram_full_size.get_proj_data_info_sptr())); + + int tested_LORs = 0; + for (int det1_downsampled = 0; det1_downsampled < 3 * 8; det1_downsampled++) + { + if (det1_downsampled % 3 == 1) + continue; // skip the central crystal of each module + for (int det2_downsampled = 0; det2_downsampled < 3 * 8; det2_downsampled++) + { + if (det2_downsampled % 3 == 1 || det1_downsampled == det2_downsampled) + continue; // skip the central crystal of each module + if (det1_downsampled / 3 == det2_downsampled / 3) + continue; // skip the LORs that lie on the same module + + int view_ds, tang_pos_ds; + pdi_downsampled->get_view_tangential_pos_num_for_det_num_pair(view_ds, tang_pos_ds, det1_downsampled, det2_downsampled); + BasicCoordinate<3, int> index_downsampled; + index_downsampled[1] = 1; // looking at central slice + index_downsampled[2] = view_ds; + index_downsampled[3] = tang_pos_ds; + + if (tang_pos_ds < pdi_downsampled->get_min_tangential_pos_num() + || tang_pos_ds > pdi_downsampled->get_max_tangential_pos_num()) + continue; + + int view_fs, tang_pos_fs; + pdi_full_size->get_view_tangential_pos_num_for_det_num_pair( + view_fs, + tang_pos_fs, + (det1_downsampled / 3) * 12 + ((det1_downsampled % 3) / 2) * 11, + (det2_downsampled / 3) * 12 + ((det2_downsampled % 3) / 2) * 11); + + BasicCoordinate<3, int> index_full_size; + index_full_size[1] = 1; // looking at central slice + index_full_size[2] = view_fs; + index_full_size[3] = tang_pos_fs; + + if (tang_pos_fs < pdi_full_size->get_min_tangential_pos_num() + || tang_pos_fs > pdi_full_size->get_max_tangential_pos_num()) + continue; + + // confirm that the difference is smaller than an empirically found value + check_if_less(std::abs(sinogram_downsampled[index_downsampled] - sinogram_full_size[index_full_size]), + 0.01, + "difference between sinogram bin is larger than expected"); + + tested_LORs++; + } + } + + info(boost::format("A total of %1% LORs were compared between the downsampled and the interpolated sinogram.") % tested_LORs); +} + void InterpolationTests::run_tests() { @@ -682,6 +959,8 @@ InterpolationTests::run_tests() scatter_interpolation_test_cyl(); scatter_interpolation_test_blocks_asymmetric(); scatter_interpolation_test_cyl_asymmetric(); + scatter_interpolation_test_blocks_downsampled(); + transaxial_upsampling_interpolation_test_blocks(); } END_NAMESPACE_STIR diff --git a/src/utilities/stir_timings.cxx b/src/utilities/stir_timings.cxx index a81fe246f3..7d65a8f4eb 100644 --- a/src/utilities/stir_timings.cxx +++ b/src/utilities/stir_timings.cxx @@ -22,6 +22,7 @@ If you want to know what is actually timed, you will have to look at the source code. */ +#include "stir/KeyParser.h" #include "stir/ProjDataInterfile.h" #include "stir/ProjDataInMemory.h" #include "stir/DiscretisedDensity.h" @@ -55,11 +56,23 @@ print_usage_and_exit() { std::cerr << "\nUsage:\nstir_timings [--name some_string] [--threads num_threads] [--runs num_runs]\\\n" << "\t[--skip-BB 1] [--skip-PP 1] [--skip-PMRT 1] [--skip-priors 1]\\\n" + << "\t[--projector_par_filename parfile]\\\n" << "\t[--image image_filename]\\\n" << "\t--template-projdata template_proj_data_filename\n\n" << "skip BB: basic building blocks; PP: Parallelproj; PMRT: ray-tracing matrix; priors: prior timing\n\n" << "Timings are reported to stdout as:\n" << "name\ttiming_name\tCPU_time_in_ms\twall-clock_time_in_ms\n"; + std::cerr << "\nExample projector-pair par-file (the following corresponds to the PMRT configuration normally used)\n" + << "projector pair parameters:=\n" + << " type := Matrix\n" + << " Projector Pair Using Matrix Parameters :=\n" + << " Matrix type := Ray Tracing\n" + << " Ray tracing matrix parameters :=\n" + << " number of rays in tangential direction to trace for each bin:= 5\n" + << " disable caching := 0\n" + << " End Ray tracing matrix parameters :=\n" + << " End Projector Pair Using Matrix Parameters :=\n" + << "End:=\n"; std::exit(EXIT_FAILURE); } @@ -89,6 +102,7 @@ class Timings : public TimedObject #ifdef STIR_WITH_Parallelproj_PROJECTOR shared_ptr parallelproj_projectors_sptr; #endif + shared_ptr parsed_projectors_sptr; shared_ptr template_proj_data_sptr; shared_ptr exam_info_sptr; shared_ptr>> objective_function_sptr; @@ -105,6 +119,7 @@ class Timings : public TimedObject } void run_it(TimedFunction f, const std::string& item, const unsigned runs = 1); + void run_projectors(const std::string& prefix, const shared_ptr proj_sptr, const unsigned runs); void run_all(const unsigned runs = 1); void init(); @@ -283,6 +298,21 @@ Timings::run_it(TimedFunction f, const std::string& item, const unsigned runs) << std::setw(24) << std::right << this->get_wall_clock_timer_value() / runs * 1000 << std::endl; } +void +Timings::run_projectors(const std::string& prefix, const shared_ptr proj_sptr, const unsigned runs) +{ + this->projectors_sptr = proj_sptr; + this->run_it(&Timings::projector_setup, prefix + "_projector_setup", 1); + this->run_it(&Timings::forward_file, prefix + "_forward_file_first", 1); + this->run_it(&Timings::forward_file, prefix + "_forward_file", runs); + this->run_it(&Timings::forward_memory, prefix + "_forward_memory", runs); + this->run_it(&Timings::back_file, prefix + "_back_file_first", 1); + this->run_it(&Timings::back_file, prefix + "_back_file", runs); + this->run_it(&Timings::back_memory, prefix + "_back_memory", runs); + this->objective_function_sptr->set_projector_pair_sptr(this->projectors_sptr); + this->run_it(&Timings::obj_func_set_up, prefix + "_LogLik set_up", 1); + this->run_it(&Timings::obj_func_grad_no_sens, prefix + "_LogLik grad_no_sens", 1); +} void Timings::run_all(const unsigned runs) { @@ -319,34 +349,18 @@ Timings::run_all(const unsigned runs) // this->objective_function.set_num_subsets(proj_data_sptr->get_num_views()/2); if (!this->skip_PMRT) { - this->projectors_sptr = this->pmrt_projectors_sptr; - this->run_it(&Timings::projector_setup, "PMRT_projector_setup", runs * 10); - this->run_it(&Timings::forward_file, "PMRT_forward_file_first", 1); - this->run_it(&Timings::forward_file, "PMRT_forward_file", 1); - this->run_it(&Timings::forward_memory, "PMRT_forward_memory", 1); - this->run_it(&Timings::back_file, "PMRT_back_file_first", 1); - this->run_it(&Timings::back_file, "PMRT_back_file", 1); - this->run_it(&Timings::back_memory, "PMRT_back_memory", 1); - this->objective_function_sptr->set_projector_pair_sptr(this->projectors_sptr); - this->run_it(&Timings::obj_func_set_up, "PMRT_LogLik set_up", 1); - this->run_it(&Timings::obj_func_grad_no_sens, "PMRT_LogLik grad_no_sens", 1); + this->run_projectors("PMRT", this->pmrt_projectors_sptr, 1); } #ifdef STIR_WITH_Parallelproj_PROJECTOR if (!skip_PP) { - this->projectors_sptr = this->parallelproj_projectors_sptr; - this->run_it(&Timings::projector_setup, "PP_projector_setup", 1); - this->run_it(&Timings::forward_file, "PP_forward_file_first", 1); - this->run_it(&Timings::forward_file, "PP_forward_file", runs); - this->run_it(&Timings::forward_memory, "PP_forward_memory", runs); - this->run_it(&Timings::back_file, "PP_back_file_first", 1); - this->run_it(&Timings::back_file, "PP_back_file", runs); - this->run_it(&Timings::back_memory, "PP_back_memory", runs); - this->objective_function_sptr->set_projector_pair_sptr(this->projectors_sptr); - this->run_it(&Timings::obj_func_set_up, "PP_LogLik set_up", 1); - this->run_it(&Timings::obj_func_grad_no_sens, "PP_LogLik grad_no_sens", 1); + this->run_projectors("PP", this->parallelproj_projectors_sptr, runs); } #endif + if (parsed_projectors_sptr) + { + this->run_projectors("parsed", this->parsed_projectors_sptr, runs); + } // write_to_file("my_timings_backproj.hv", *this->image_sptr); if (!skip_priors) @@ -445,6 +459,7 @@ main(int argc, char** argv) std::string image_filename; std::string template_proj_data_filename; + std::string projector_par_filename; std::string prog_name = argv[0]; unsigned num_runs = 3; int num_threads = get_default_num_threads(); @@ -477,6 +492,8 @@ main(int argc, char** argv) skip_PP = std::atoi(argv[1]) != 0; else if (!strcmp(argv[0], "--skip-priors")) skip_priors = std::atoi(argv[1]) != 0; + else if (!strcmp(argv[0], "--projector_par_filename")) + projector_par_filename = argv[1]; else print_usage_and_exit(); argv += 2; @@ -495,6 +512,15 @@ main(int argc, char** argv) timings.skip_PMRT = skip_PMRT; timings.skip_PP = skip_PP; timings.skip_priors = skip_priors; + if (!projector_par_filename.empty()) + { + KeyParser parser; + parser.add_start_key("Projector pair parameters"); + parser.add_parsing_key("type", &timings.parsed_projectors_sptr); + parser.add_stop_key("END"); + if (!parser.parse(projector_par_filename.c_str())) + error("Error parsing " + projector_par_filename); + } timings.run_all(num_runs); return EXIT_SUCCESS;