diff --git a/.appveyor.yml b/.appveyor.yml index 6125a65d9..326e4515a 100644 --- a/.appveyor.yml +++ b/.appveyor.yml @@ -3,6 +3,7 @@ shallow_clone: true skip_commits: files: - .github/**/* + - 'CITATION.cff' - '**/*.md' - '**/*.html' - '**/*.htm' diff --git a/.gitattributes b/.gitattributes index 7d845e399..dc8338a6a 100644 --- a/.gitattributes +++ b/.gitattributes @@ -12,6 +12,8 @@ *.v -text -diff *.s -text -diff *.scn -text -diff +*.l -text -diff +*.root -text -diff *.gz -text -diff *.nii -text -diff *.sh eol=lf diff --git a/.github/PULL_REQUEST_TEMPLATE.md b/.github/PULL_REQUEST_TEMPLATE.md new file mode 100644 index 000000000..aed2641f0 --- /dev/null +++ b/.github/PULL_REQUEST_TEMPLATE.md @@ -0,0 +1,22 @@ + + +## Changes in this pull request + + +## Testing performed + + +## Related issues + + + +## Checklist before requesting a review + + - [] I have performed a self-review of my code + - [] I have added docstrings/doxygen in line with the guidance in the developer guide + - [] I have implemented unit tests that cover any new or modified functionality (if applicable) + - [] The code builds and runs on my machine + - [] `documentation/release_XXX.md` has been updated with any functionality change (if applicable) diff --git a/.github/workflows/build-test.yml b/.github/workflows/build-test.yml index 128528d99..1be05e812 100644 --- a/.github/workflows/build-test.yml +++ b/.github/workflows/build-test.yml @@ -7,6 +7,7 @@ on: - tof_sino_UCL paths-ignore: - '.appveyor.yml' + - 'CITATION.cff' - '**/*.md' - '**/*.html' - '**/*.htm' @@ -18,6 +19,7 @@ on: - tof_sino_UCL paths-ignore: - '.appveyor.yml' + - 'CITATION.cff' - '**/*.md' - '**/*.html' - '**/*.htm' @@ -334,7 +336,7 @@ jobs: ./run_scatter_tests.sh ./run_test_zoom_image.sh ./run_ML_norm_tests.sh - if [[ $BUILD_FLAGS == *"DDISABLE_CERN_ROOT=0"* ]]; then ./run_root_GATE.sh; fi + if test "${{matrix.ROOT}}XX" == "ONXX"; then ./run_root_GATE.sh; fi ./run_tests_modelling.sh cd ${GITHUB_WORKSPACE}/recon_test_pack/SPECT ./run_SPECT_tests.sh @@ -345,14 +347,17 @@ jobs: if: failure() with: name: recon_test_pack_log_files-${{ matrix.os }}-${{ matrix.compiler }}${{ matrix.compiler_version }}-${{ matrix.BUILD_TYPE }}-pp=${{ matrix.parallelproj }}-ROOT=${{ matrix.ROOT }} - path: ${{ github.workspace }}/recon_test_pack/**/*.log + path: | + ${{ github.workspace }}/recon_test_pack/**/*.log + ${{ github.workspace }}/recon_test_pack/**/my_*v + ${{ github.workspace }}/recon_test_pack/**/my_*s retention-days: 7 # 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' && github.event.inputs.debug_enabled }} + if: ${{ github.event_name == 'workflow_dispatch' && inputs.debug_enabled == 'true' }} - name: examples shell: bash diff --git a/.github/workflows/cffconvert.yml b/.github/workflows/cffconvert.yml new file mode 100644 index 000000000..707a71c4b --- /dev/null +++ b/.github/workflows/cffconvert.yml @@ -0,0 +1,19 @@ +name: cffconvert + +on: + push: + paths: + - CITATION.cff + +jobs: + validate: + name: "validate" + runs-on: ubuntu-latest + steps: + - name: Check out a copy of the repository + uses: actions/checkout@v2 + + - name: Check whether the citation metadata from CITATION.cff is valid + uses: citation-file-format/cffconvert-github-action@2.0.0 + with: + args: "--validate" diff --git a/.mailmap b/.mailmap index e0b863e5c..66b9eb55f 100644 --- a/.mailmap +++ b/.mailmap @@ -45,8 +45,10 @@ Palak Wadhwa Rebecca Gillen <43674917+RebeccaGillen@users.noreply.github.com> Richard Brown <33289025+rijobro@users.noreply.github.com> Richard Brown -Robert Twyman -Robert Twyman +Robert Twyman Skelly +Robert Twyman Skelly +Robert Twyman Skelly +Robert Twyman Skelly <117300855+Robert-PrescientImaging@users.noreply.github.com> Yu-jung Tsai Gemma Fardell Gemma Fardell <47746591+gfardell@users.noreply.github.com> @@ -62,5 +64,6 @@ Nicole Jurjew Tahereh Niknejad Tahereh Niknejad Sam D Porter <92305641+samdporter@users.noreply.github.com> +Sam D Porter <92305641+samdporter@users.noreply.github.com> Matthew Strugari Matthew Strugari <56315593+mastergari@users.noreply.github.com> diff --git a/.zenodo.json b/.zenodo.json deleted file mode 100644 index ce2659bc9..000000000 --- a/.zenodo.json +++ /dev/null @@ -1,89 +0,0 @@ -{ - "title": "STIR: Software for Tomographic Image Reconstruction", - "keywords": [ - "image-reconstruction", "pet", "spect", "medical-imaging", "open-source software"], - "related_identifiers": [ - {"identifier": "10.1088/0031-9155/57/4/867", "relation": "documents"}, - {"identifier": "10.1186/s40658-019-0248-9", "relation": "documents"}, - {"identifier": "10.3390/jimaging8060172", "relation": "documents"}, - {"identifier": "10.1109/NSSMIC.2018.8824341", "relation": "documents"}, - {"identifier": "10.1186/2197-7364-1-s1-a44", "relation": "documents"}, - {"identifier": "10.1109/nssmic.2013.6829258", "relation": "documents"}, - {"identifier": "10.1118/1.4816676", "relation": "documents"}, - {"identifier": "10.1007/s12149-011-0514-y", "relation": "documents"}, - {"identifier": "10.1016/j.compmedimag.2011.01.002", "relation": "documents"}, - {"identifier": "10.1088/1742-6596/317/1/012002", "relation": "documents"}, - {"identifier": "10.1109/nssmic.2008.4774198", "relation": "documents"}, - {"identifier": "10.1109/nssmic.2006.354345", "relation": "documents"}, - {"identifier": "10.1109/nssmic.2005.1596753", "relation": "documents"}, - {"identifier": "10.1109/nssmic.2004.1466455", "relation": "documents"}, - {"identifier": "10.1109/nssmic.2004.1466782", "relation": "documents"}, - {"identifier": "10.1109/nssmic.2002.1239610", "relation": "documents"}, - {"identifier": "10.1109/NSSMIC.2001.1008688", "relation": "documents"}, - {"identifier": "10.1088/0031-9155/45/8/325", "relation": "documents"}, - {"identifier": "10.1007/978-3-642-60125-5_50", "relation": "documents"}, - {"identifier": "10.1088/1361-6420/ab013f", "relation": "documents"}, - {"identifier": "10.1109/TRPMS.2018.2884176", "relation": "documents"}, - {"identifier": "10.1109/nssmic.2018.8824312", "relation": "documents"}, - {"identifier": "10.1109/NSS/MIC42101.2019.9059665", "relation": "documents"}], - "communities": [{"identifier": "ccp-petmr"},{"identifier": "synerbi"}], - "creators": [ - {"name": "Mustafovic, Sanida", "affiliation": "Imperal College London (UK)"}, - {"name": "Efthimiou, Nikos", "orcid": "0000-0003-1947-5033"}, - {"name": "Brown, Richard", "orcid": "0000-0001-6989-9200", "affiliation": "University College London"}, - {"name": "Tsoumpas, Charalampos"}, - {"name": "Twyman, Robert", "affiliation": "University College London"}, - {"name": "Deidda, Daniel", "orcid": "0000-0002-2766-4339"}, - {"name": "Borgeaud, Tim", "affiliation": "Hammersmith Imanet Ltd"}, - {"name": "Falcon, Carles"}, - {"name": "Khateri, Parisa", "affiliation": "ETH Zuerich"}, - {"name": "Wadhwa, Palak", "affiliation": "University of Leeds (UK)"}, - {"name": "Beisel, Tobias"}, - {"name": "Strugari, Matthew", "affiliation": "Dalhousie University (Canada)"}, - {"name": "Jacobson, Matthew"}, - {"name": "Gillman, Ashley", "orcid": "0000-0001-9130-1092"}, - {"name": "Zverovich, Alexey", "affiliation": "Brunel University (UK)"}, - {"name": "Emond, Elise", "affiliation": "University College London"}, - {"name": "Biguri, Ander", "orcid": "0000-0002-2636-3032"}, - {"name": "Fuster Marti, Berta", "affiliation": "University of Barcelona (Spain)"}, - {"name": "Labbe, Claire"}, - {"name": "Fischer, Jannis", "orcid": "0000-0002-8329-0220"}, - {"name": "Jehl, Markus", "affiliation":"Positrigo"}, - {"name": "Roethlisberger, Michael", "affiliation": "ETH Zuerich"}, - {"name": "Aguiar, Pablo"}, - {"name": "Brusaferri, Ludovica", "affiliation": "University College London"}, - {"name": "Bertolli, Ottavia", "affiliation": "University College London"}, - {"name": "Pasca, Edoardo", "orcid": "0000-0001-6957-2160"}, - {"name": "Niknejad, Tahereh"}, - {"name": "Dikaios, Nikos"}, - {"name": "Sadki, Mustapha", "affiliation": "Brunel University (UK)"}, - {"name": "Fardell, Gemma", "orcid": "0000-0003-2388-5211", "affiliation":"UK Research & Innovation"}, - {"name": "Kerrouche, Nacer", "affiliation": "Hammersmith Imanet Ltd"}, - {"name": "Ovtchinnikov, Evgueni"}, - {"name": "Ehrhardt, Matthias J.", "orcid": "0000-0001-8523-353X"}, - {"name": "Schmidtlein, C. Ross"}, - {"name": "Valente, Patrick", "affiliation": "Brunel University (UK)"}, - {"name": "Thomas, Benjamin", "orcid": "0000-0002-9784-1177"}, - {"name": "Schramm, Georg", "affiliation":"Katholieke Universiteit Leuven (Belgium)"}, - {"name": "Völgyes, David"}, - {"name": "Dinelle, Katie"}, - {"name": "Belluzzo, Damiano", "affiliation": "Hospedale San Raffaele Milano (Italy)"}, - {"name": "Ching, Daniel"}, - {"name": "Hague, Darren", "affiliation": "Brunel University (UK)"}, - {"name": "Tunnicliffe, Harry", "affiliation":"University of Leeds (UK)"}, - {"name": "Chen, Gefei"}, - {"name": "Dao, Viet Anh", "affiliation":"University of Leeds (UK)"}, - {"name": "Mikhaylova, Ekaterina", "affiliation":"Positrigo"}, - {"name": "da Costa-Luis, Casper O.", "orcid": "0000-0002-7211-1557"}, - {"name": "Porter, Sam D.", "affiliation": "University College London"}, - {"name": "Rashidnasab, Alaleh", "affiliation": "University College London"}, - {"name": "Whitehead, Alexander C.", "affiliation": "University College London"}, - {"name": "Gillen, Rebecca", "affiliation": "University College London"}, - {"name": "Vavrek, Jayson"} - {"name": "Kohr, Holger"}, - {"name": "Tsai, Yu-jung", "affiliation": "University College London"}, - {"name": "Kohr, Holger"}, - {"name": "tokkot"}, - {"name": "El Katib, Mahmoud"}, - {"name": "Thielemans, Kris", "orcid": "0000-0002-5514-199X", "affiliation": "University College London"}] -} diff --git a/CITATION.cff b/CITATION.cff new file mode 100644 index 000000000..77d74b039 --- /dev/null +++ b/CITATION.cff @@ -0,0 +1,229 @@ +# This CITATION.cff file was generated with cffinit. +# Visit https://bit.ly/cffinit to generate yours today! + +cff-version: 1.2.0 +title: STIR Software for Tomographic Image Reconstruction +message: >- + If you use this software, please cite it using the + metadata from this file. +type: software +authors: + - family-names: Mustafovic + given-names: Sanida + affiliation: Imperial College London (UK) + - family-names: Efthimiou + given-names: Nikos + orcid: 'https://orcid.org/0000-0003-1947-5033' + - family-names: Brown + given-names: Richard + orcid: 'https://orcid.org/0000-0001-6989-9200' + affiliation: University College London + - family-names: Twyman Skelly + given-names: Robert + affiliation: University College London and Prescient Imaging + - family-names: Tsoumpas + given-names: Charalampos + - family-names: Deidda + given-names: Daniel + orcid: 'https://orcid.org/0000-0002-2766-4339' + affiliation: National Physics Laboratory (UK) + - family-names: Falcon + given-names: Carles + - family-names: Borgeaud + given-names: Tim + affiliation: Hammersmith Imanet Ltd + - family-names: Khateri + given-names: Parisa + affiliation: ETH Zuerich + - family-names: Jehl + given-names: Markus + affiliation: Positrigo + - family-names: Wadhwa + given-names: Palak + affiliation: University of Leeds (UK) + - family-names: Strugari + given-names: Matthew + affiliation: Dalhousie University (Canada) + - family-names: Beisel + given-names: Tobias + - family-names: Jacobson + given-names: Matthew + - family-names: Gillman + given-names: Ashley + orcid: 'https://orcid.org/0000-0001-9130-1092' + affiliation: >- + Commonwealth Scientific and Industrial Research + Organisation, and University of Queensland + - family-names: Zverovich + given-names: Alexey + affiliation: Brunel University (UK) + - family-names: Emond + given-names: Elise + affiliation: University College London + - family-names: Biguri + given-names: Ander + orcid: 'https://orcid.org/0000-0002-2636-3032' + affiliation: University College London + - family-names: Fuster Marti + given-names: Berta + affiliation: University of Barcelona (Spain) + - family-names: Labbe + given-names: Claire + - family-names: Fischer + given-names: Jannis + orcid: 'https://orcid.org/0000-0002-8329-0220' + - family-names: Roethlisberger + given-names: Michael + affiliation: ETH Zuerich + - family-names: Aguiar + given-names: Pablo + - family-names: Brusaferri + given-names: Ludovica + affiliation: University College London + - family-names: Pasca + given-names: Edoardo + affiliation: UK Research & Innovation + orcid: 'https://orcid.org/0000-0001-6957-2160' + - family-names: Thomas + given-names: Benjamin + orcid: 'https://orcid.org/0000-0002-9784-1177' + affiliation: University College London + - family-names: Bertolli + given-names: Ottavia + affiliation: University College London + - family-names: Niknejad + given-names: Tahereh + - family-names: Dikaios + given-names: Nikos + - family-names: Sadki + given-names: Mustapha + affiliation: Brunel University (UK) + - family-names: Fardell + given-names: Gemma + orcid: 'https://orcid.org/0000-0003-2388-5211' + affiliation: UK Research & Innovation + - family-names: Kerrouche + given-names: Nacer + affiliation: Hammersmith Imanet Ltd + - family-names: Ovtchinnikov + given-names: Evgueni + orcid: 'https://orcid.org/0000-0002-9359-6514' + affiliation: UK Research & Innovation + - family-names: Ehrhardt + given-names: Matthias J. + orcid: 'https://orcid.org/0000-0001-8523-353X' + - family-names: Schmidtlein + given-names: C Ross + - family-names: Valente + given-names: Patrick + affiliation: Brunel University (UK) + - family-names: Schramm + given-names: Georg + affiliation: Katholieke Universiteit Leuven (Belgium) + - family-names: Völgyes + given-names: David + - family-names: Dinelle + given-names: Katie + - family-names: Belluzzo + given-names: Damiano + affiliation: Hospedale San Raffaele Milano (Italy) + - family-names: Ching + given-names: Daniel + - family-names: Hague + given-names: Darren + affiliation: Brunel University (UK) + - family-names: Tunnicliffe + given-names: Harry + affiliation: University of Leeds (UK) + - family-names: Chen + given-names: Gefei + - family-names: Dao + given-names: Viet Anh + affiliation: University of Leeds (UK) + - family-names: Mikhaylova + given-names: Ekaterina + affiliation: Positrigo + - family-names: Porter + given-names: Sam David + affiliation: University College London, National Physics Laboratory (UK) + - family-names: da Costa-Luis + given-names: Casper O. + orcid: 'https://orcid.org/0000-0002-7211-1557' + - family-names: Rashidnasab + given-names: Alaleh + affiliation: University College London + - family-names: Whitehead + given-names: Alexander C. + affiliation: University College London + - family-names: Gillen + given-names: Rebecca + affiliation: University College London + - family-names: Vavrek + given-names: Jayson + - family-names: Tsai + given-names: Yu-jung + affiliation: University College London + - family-names: Kohr + given-names: Holger + - name: tokkot + - family-names: El Katib + given-names: Mahmoud + - family-names: Thielemans + given-names: Kris + orcid: 'https://orcid.org/0000-0002-5514-199X' + affiliation: University College London +identifiers: + - type: doi + value: 10.1088/0031-9155/57/4/867 + - type: doi + value: 10.1186/s40658-019-0248-9 + - type: doi + value: 10.3390/jimaging8060172 + - type: doi + value: 10.1109/NSSMIC.2018.8824341 + - type: doi + value: 10.1186/2197-7364-1-s1-a44 + - type: doi + value: 10.1109/nssmic.2013.6829258 + - type: doi + value: 10.1118/1.4816676 + - type: doi + value: 10.1007/s12149-011-0514-y + - type: doi + value: 10.1016/j.compmedimag.2011.01.002 + - type: doi + value: 10.1088/1742-6596/317/1/012002 + - type: doi + value: 10.1109/nssmic.2008.4774198 + - type: doi + value: 10.1109/nssmic.2006.354345 + - type: doi + value: 10.1109/nssmic.2005.1596753 + - type: doi + value: 10.1109/nssmic.2004.1466455 + - type: doi + value: 10.1109/nssmic.2004.1466782 + - type: doi + value: 10.1109/nssmic.2002.1239610 + - type: doi + value: 10.1109/NSSMIC.2001.1008688 + - type: doi + value: 10.1088/0031-9155/45/8/325 + - type: doi + value: 10.1007/978-3-642-60125-5_50 + - type: doi + value: 10.1088/1361-6420/ab013f + - type: doi + value: 10.1109/TRPMS.2018.2884176 + - type: doi + value: 10.1109/nssmic.2018.8824312 + - type: doi + value: 10.1109/NSS/MIC42101.2019.9059665 +repository-code: 'https://github.com/UCL/STIR' +keywords: + - pet + - spect + - medical imaging + - image reconstruction + - open-source software +license: Apache-2.0 diff --git a/CMakeLists.txt b/CMakeLists.txt index cec25b64b..fc7b94e51 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -55,10 +55,10 @@ else() endif() ####### Set Version number etc -set(VERSION_MAJOR 5) -set(VERSION_MINOR 1) -set(VERSION_PATCH 2) -set(VERSION 050102) # only used in STIRConfig.h.in +set(VERSION_MAJOR 6) +set(VERSION_MINOR 0) +set(VERSION_PATCH 0) +set(VERSION 060000) # only used in STIRConfig.h.in and swig/CMakeLists.txt set(STIR_VERSION ${VERSION_MAJOR}.${VERSION_MINOR}.${VERSION_PATCH}) diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 7210cb764..3c8ea7411 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -44,7 +44,6 @@ Please by mindful about the resources used by our Continuous Integration (CI) wo - Use specific keywords in the first line of the last commit that you push to prevent CI being run: - `[ci skip]` skips all CI runs (e.g. when you only change documentation, or when your update isn't ready yet) - `[actions skip]` does not run GitHub Actions, see [here](https://github.blog/changelog/2021-02-08-github-actions-skip-pull-request-and-push-workflows-with-skip-ci/). Note: this can be in the main commit message. - - `[travis skip]` does not run Travis-CI, see [here](https://docs.travis-ci.com/user/customizing-the-build/#skipping-a-build). Note: this can be in the main commit message. - `[skip appveyor]` does not run Appveyor, see [here](https://www.appveyor.com/docs/how-to/filtering-commits/#skip-directive-in-commit-message) 8. After acceptance of your PR, go home with a nice warm feeling. diff --git a/README.md b/README.md index c5213623a..65096baa2 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,10 @@ # STIR: Software for Tomographic Image Reconstruction. +[![GitHub Actions status](https://github.com/UCL/STIR/actions/workflows/build-test.yml/badge.svg)](https://github.com/UCL/STIR/actions/workflows/build-test.yml) +[![Appveyor Build status](https://ci.appveyor.com/api/projects/status/ga9xd1vsy0ik1soq/branch/master?svg=true)](https://ci.appveyor.com/project/KrisThielemans/stir/branch/master) +[![Codacy Badge](https://app.codacy.com/project/badge/Grade/c561421aa1af41e4a9ef779003f6fff0)](https://app.codacy.com/gh/UCL/STIR/dashboard?utm_source=gh&utm_medium=referral&utm_content=&utm_campaign=Badge_grade) +[![Zenodo DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.6604468.svg)](https://doi.org/10.5281/zenodo.6604468) + STIR is Open Source software for use in tomographic imaging. Its aim is to provide a Multi-Platform Object-Oriented framework for all data manipulations in tomographic imaging. Currently, the emphasis is on @@ -17,10 +22,3 @@ information. This software is distributed under an open source license, see [LICENSE.txt](LICENSE.txt) for details. -## Build and test status of the master branch -- Travis (tests Linux and MacOS) -[![Travis Build Status](https://travis-ci.org/UCL/STIR.svg?branch=master)](https://travis-ci.org/UCL/STIR) - -- Appveyor (tests Windows) -[![Appveyor Build status](https://ci.appveyor.com/api/projects/status/ga9xd1vsy0ik1soq/branch/master?svg=true)](https://ci.appveyor.com/project/KrisThielemans/stir/branch/master) - diff --git a/VERSION.txt b/VERSION.txt index 61fcc8735..91ff57278 100644 --- a/VERSION.txt +++ b/VERSION.txt @@ -1 +1 @@ -5.1.2 +5.2.0 diff --git a/documentation/STIR-UsersGuide.tex b/documentation/STIR-UsersGuide.tex index 80a91a9ae..0d14db36b 100644 --- a/documentation/STIR-UsersGuide.tex +++ b/documentation/STIR-UsersGuide.tex @@ -44,7 +44,7 @@ \\[3cm] \textbf{{\huge User's Guide\\ - Version 5.1}} + Version 5.2}} \end{center} \end{spacing} @@ -53,7 +53,7 @@ \noindent K. Thielemans \\ -{\it \small Hammersmith Imanet Ltd; Algorithms and Software Consulting Ltd; University College London}\\ +{\it \small University College London; Algorithms and Software Consulting Ltd (formerly Hammersmith Imanet Ltd}\\ Ch. Tsoumpas\\ {\it \small Hammersmith Imanet Ltd; Imperial College London; King's College London}\\ D. Sauge, C. Labb\'e, C. Morel \\ @@ -142,11 +142,42 @@ \section{ Any comments on documentation, and especially contributions are always welcome. Please use the stir-users mailing list. +\section{Installing STIR via a pre-built version } +\subsection{ +Installing via conda} +The easiest way is to use \texttt{conda} (or one of its derivatives such as \texttt{mamba}): +\begin{verbatim} +conda config --add channels conda-forge +conda create -n stirenv stir matplotlib +conda activate stirenv + +# use STIR executable or python + +# deactivate the environment +conda deactivate +\end{verbatim} +See \url{https://stir.sourceforge.net/wiki/index.php/Installing_STIR_with_conda}{our Wiki page} for +more detail, including on how to get a development version. + +\subsection{SIRF distributions} +CCP SyneRBI distributes SIRF with its dependencies, including STIR, +see \url{https://www.ccpsynerbi.ac.uk/downloads/}{its download page} for +a Virtual Machine, docker etc. However, only the STIR libraries are included in +these distributions. You can easily install ``full'' STIR on the Virtual Machine though, +see \url{https://github.com/SyneRBI/SIRF-SuperBuild/blob/master/VirtualBox/documentation/README.md}{its documentation}. + \section{ -Installation} +Installation from source} This section describes how to install STIR. It is complemented by information on the \url{http://stir.sourceforge.net/wiki}{STIR Wiki} with information for specific systems etc. +Note that as opposed to using the instructions below, you can also +the \url{https://github.com/SyneRBI/SIRF-SuperBuild/blob/master/README.md}{SIRF-Superbuild} +of \url{https://www.ccpsynerbi.ac.uk/}{SyneRBI} to build STIR with CMake. This will take care of all dependencies. +Note that you will need to use advanced CMake variables \texttt{STIR\_BUILD\_EXECUTABLES=ON} +and \texttt{STIR\_BUILD\_SWIG\_PYTHON=ON}. + + \subsection{ Installing source files} @@ -156,26 +187,15 @@ \subsection{ \end{center} or from \url{https://github.com/UCL/STIR/releases}{github.com/UCL/STIR/releases}.\\ Alternatively, if you are feeling adventurous, you can get the most-recent -developer version (no guarantees!) of STIR from\\ +developer version (no guarantees!) of STIR from \url{https://github.com/UCL/STIR}{https://github.com/UCL/STIR}.\\ -Download the source files in the zip format. You can then use unzip from \\ -\url{http://www.info-zip.org/pub/infozip/}{http://www.info-zip.org/pub/infozip/}. -Extract using -\cmdline{unzip -a file.zip} - - -The \texttt{-a} option makes sure that files are extracted using the -end-of-line convention appropriate for your system. Note that -other programs that can handle zip files (such as WinZip) should -work as well, although you might have problems with the EOL convention. - Note that you can put the distribution in any location (it does not have to be your home directory). -The result is a \textbf{STIR} directory and subdirectories, described -in the annex section. +After unpacking, you should have a \textbf{STIR} directory (subdirectories are described +in the doxygen documentation). \subsection{ Installing external software} @@ -183,6 +203,13 @@ \subsection{ On most systems, you should be able to get these using your package manager. Please check the Wiki for most up-to-date information. +\subsubsection{ +C++ Compiler} + +In order to compile and run the \textit{STIR} programs, you will need a compiler +and associated tools such as \texttt{CMake}. These days, any compiler should work. We would love +to hear from any attempts where there are failures. + \subsubsection{BOOST} The only required external library is the well-respected \textit{boost} library. If it is not installed on your system, you can download it from @@ -195,29 +222,19 @@ \subsubsection{BOOST} \subsubsection{JSON} The \texttt{HUToMu} and \texttt{Radionuclide} database functionality parse a file in JSON format. We use the \url{https://github.com/nlohmann/json}{nlohman\_json} library for this. Installation -instructions for this library are available on its github page. +instructions for this library are available on its github page or you can get it via \texttt{conda}. +Without it, we only support a few basic radionuclides. \subsubsection{ -C++ Compiler} +Enabling ECAT 7 support\label{sec:ECAT67support}} +{\em ECAT7 support is no longer tested.} -In order to compile and run the \textit{STIR} programs, you will need a compiler -and associated tools. These days, any compiler should work. See the -\url{http://sourceforge.net/stir}{STIR Wiki} for more -specific information of which compilers we have tried. -We would love -to hear from any attempts of using another compiler. - -\subsubsection{ -Enabling ECAT 7 support} -\label{sec:ECAT67support} Older CTI (Siemens) scanners use a file format called ECAT\texttrademark{}. At present, \textit{STIR} uses parts of the Louvain la Neuve ecat library (called LLN in the rest of this document).\footnote{\textit{STIR} versions 1.? came with specific files for ECAT6 support without the need for the LLN library. However, due to license restrictions this is now no longer the case.} -The library might still be available on -\url{ftp://ftp.topo.ucl.ac.be/pub/ecat }{ftp://ftp.topo.ucl.ac.be/pub/ecat}. -It is also available for GATE users. +The library could be available on the OpenGATE website. You have to download that library and issue 'make --f Makefile.unix' (or 'make --f Makefile.cygwin' if you are using CYGWIN on Windows) @@ -235,10 +252,26 @@ \subsubsection{ For most operating systems this can be done via your package manager which we highly recommend. You could also download from the \url{https://www.hdfgroup.org/downloads/hdf5/}{HDF5 group download page}. +\subsubsection{ +Enabling ROOT support\label{sec:installROOT}} +STIR can read CERN ROOT files from GATE (see section \ref{sec:ROOTIO}). Check the \url{https://root.cern/install/}{installation instructions for ROOT}. + +\subsubsection{ +Enabling ITK support\label{sec:installITK}} +STIR can use \url{http://www.itk.org}{ITK}, a large open source library. Specifying this + enables NRRD, MetaIO and Nifti IO, and DICOM reading. + Use your package manager, including conda or pip, or check the \url{https://itk.org/download/}{installation instructions for ITK}. + +\subsubsection{ +Enabling AVW support\label{sec:installAVW}} +AVW was the library used by the commercial program \url{https://analyzedirect.com/}{Analyze}. +\textit{AVW support has not been tested since about 2010 and we have no idea of STIR still builds with it.} AVW support will be dropped in STIR 6.0. +See \ref{sec:convAVW} for usage. + \subsection{ Building} Since version 2.2, \textit{STIR} contains files to build STIR using the platform-independent -\url{"http://www.cmake.org"}{CMake}. This is now the only to build STIR as +\url{"http://www.cmake.org"}{CMake}. This is now the only option to build STIR as it is easier for configuration and finding system dependencies. \subsubsection{ Using CMake} @@ -309,10 +342,6 @@ \subsubsection{ % ROW {\raggedright \textit{STIR\_OPENMP}} & {\raggedright Toggles between ON, OFF. Enable threaded processing using OPENMP. -\footnote{Since version 2.0, \textit{STIR} contains preliminary code for running \texttt{FBP2D} in threaded mode -(contribution by Tobias Beisel, Univ. of Paderborn). You need a compiler that -supports \texttt{OPENMP}. Note however, that this code doesn't presently result -in a decent speed-up of \texttt{FBP2D}.} See section \ref{sec:RunningWithOPENMP} for information on how to run programs using OPENMP. } \\ \hline @@ -347,14 +376,13 @@ \subsubsection{ capabilities of STIR. \begin{itemize} \item LLN files for ECAT support via \texttt{LLN\_INCLUDE\_DIRS} and \texttt{LLN\_LIBRARIES}. -\item AVW\texttrademark{} (a commercial library\footnote{See \url{http://www.mayo.edu/bir/Software/AVW/AVW1.html} -{www.mayo.edu/bir/Software/AVW/AVW1.html}.}) -via \texttt{AVW\_ROOT\_DIR}. See section \ref{sec:convAVW} for usage.\\ -\textit{AVW support has not been tested since about 2010.} -\item ITK\_DIR, use IO from \url{http://www.itk.org}{ITK}, a large open source library. Specifying this - enables NRRD, MetaIO and Nifti IO. +See section \ref{sec:ECAT67support}. +\item AVW\texttrademark{} [deprecated] (a commercial library, probably no longer available) +via \texttt{AVW\_ROOT\_DIR}. See section \ref{sec:installAVW}. +\item ITK\_DIR, use IO from ITK, see section \ref{sec:installITK}. \item CERN ROOT files used by GEANT and GATE via \texttt{ROOT\_DIR}. If this is not set, we will try to get it from the \texttt{ROOTSYS} environment (or CMake) variable. + See \ref{sec:installROOT}. \item GE RDF\texttrademark{} 9 support via \texttt{HDF5\_ROOT} and related variables (requires the HDF5 library). \end{itemize} @@ -371,14 +399,15 @@ \subsubsection{ \item \textit{CMAKE\_CXX\_STANDARD} can be used to tell CMake to add flags to your compiler to use a particular version of the C++ standard (if it supports it). It is set to \texttt{11} since STIR 5.0. Other possible allowed values are -\texttt{14} and \texttt{17}. More recent version have not yet been tested. +\texttt{14}, \texttt{17} and \textit{20}. More recent version have not yet been tested. \end{itemize} When building the Python code, and you have multiple versions of Python installed, it is often necessary to specify the correct Python executable that you want, together with -the libraries. You can use the \texttt{PYTHON\_EXECUTABLE} variable for this. If you -specify the full path to the actual executable, this should be sufficient. If not, set -\texttt{PYTHON\_LIBRARY} as well. +the libraries. You can use the \texttt{Python\_EXECUTABLE} variable for this. If you +specify the full path to the actual executable, this should be sufficient. +\footnote{If you are using CMake 3.13 or older, use \texttt{PYTHON\_EXECUTABLE}. +In this case, you might have to set \texttt{PYTHON\_LIBRARY} as well.} There are various other variables that can be set, some of which are only visible if you toggled the display of \texttt{Advanced} variables on. They should only be set by advanced users, or if a package @@ -5723,12 +5752,12 @@ \section{Using STIR in an external C++ project \label{sec:ExternalProjectC++}} STIR exports its CMake settings. Therefore, an external project can do \begin{verbatim} -find_package(STIR 5.1 CONFIG) +find_package(STIR 5.2 CONFIG) add_library(my_lib file1.cxx file2.cxx) # my_lib uses STIR functionality -target_link_libraries(my_lib ${STIR_LIBRARIES}) +target_link_libraries(my_lib PUBLIC ${STIR_LIBRARIES}) add_executable(my_exe my_exe.cxx ${STIR_REGISTRIES}) -target_link_libraries(my_exe my_lib) +target_link_libraries(my_exe PUBLIC my_lib) \end{verbatim} In addition, if your CMake is older than 3.12, you need to add \begin{verbatim} @@ -5749,20 +5778,9 @@ \section{ The \textit{STIR} library, in its current state, possesses many capabilities. The developers, however, look forward to still further increases -in the flexibility and power of the software. Some of the developments -being discussed are: - -\begin{itemize} -\item -expanded library of polymorphic classes (e.g. image grids, and -ordered subsets) -\item -additional scanners, also for SPECT -\item -additional data formats support, without conversion to Interfile. -\item -point-spread function reconstruction -\end{itemize} +in the flexibility and power of the software. +Please check the \url{https://github.com/UCL/STIR/milestones}{GitHub milestones} +for some information. While support for the library is on a voluntary basis, users @@ -5771,35 +5789,12 @@ \section{ \url{http://stir.sourceforge.net/ }{http://stir.sourceforge.net}) where they can follow developments of the software and obtain helpful information from other users. Questions will ONLY be -answered (if at all) when directed to the mailing list. +answered (if at all) when creating a \url{https://github.com/UCL/STIR/issues}{GitHub issue} (preferred) or +when directed to the mailing list. Commercial support is available from \url{http://asc.uk.com}{Algorithms and Software Consulting Ltd}. -Below, we list of some of the features that might make it into the next releases. -However, which feature is actually finalised/implemented depends -on the needs of the developers. If you want one of these features -and are willing to help, let us know. -\begin{itemize} -\item -More automatic testing programs -\item -More algorithms: potentially ART, -OSCB [Ben99b] -\item -More projectors -\item -More priors -\item -Extending the parallelisation of OSMAPOSL and OSSPS to FBP3DRP etc, or using OPEMMP -\item -Compatibility of the interpolating backprojector with recent -processors. -\item -More kinetic models: Spectral Analysis, Logan Plot -\end{itemize} - - diff --git a/documentation/STIR-developers-overview.tex b/documentation/STIR-developers-overview.tex index b06805c05..35a6d4370 100644 --- a/documentation/STIR-developers-overview.tex +++ b/documentation/STIR-developers-overview.tex @@ -16,7 +16,7 @@ \textbf{{\huge STIR \\ Overview for developers}}\\ \textbf{Kris Thielemans}\\ -\textbf{\textit{version 5.1}} +\textbf{\textit{version 5.2}} \end{center} @@ -60,7 +60,8 @@ \section{ doxygen on the source files that come in the STIR distribution). \subsection{Language support} -STIR is written in C++ and currently requires C++-11, but it should be compatible with nNewer versions of the C++ standard. +STIR is written in C++ and currently requires C++-11, but it is compatible with newer versions of the C++ standard. +We will enforce C++-14 from STIR 6.0. Python and MATLAB support is provided via \R2Lurl{http://www.swig.org/}{SWIG}. This means that Python/MATLAB interfaces follow the C++ classes closely, although some differences are required as these languages do not support templates for instance. @@ -139,31 +140,25 @@ \subsection{ In 3D PET, two data storage modes are generally used for a \textbf{segment}\footnote{{The -GE Advance file format does \textit{not} store the data per segment, -but per view. In addition, the segments are then mixed. However, +GE file formats does \textit{not} store the data per segment, +but per view. However, once read into STIR, this organisation is no longer available. -\\ -The ECAT6 sinograms also have no easy way to understand which -ring-pair corresponds to which sinogram number for 3D PET. This -is the main reason why STIR does not attempt to read ECAT6 sinograms -directly, but leaves this to a separate conversion utility.}} +}} \begin{itemize} \item where the 3D data is ordered by \textbf{sinogram} (i.e \textbf{axial position}, \textbf{view -angle}, \textbf{tangential position}) (CTI terminology: \textit{Volume} mode) +angle}, \textbf{tangential position}) (see Figure \ref{fig:sinogramstorage}) \item where the 3D dataset is ordered by \textbf{view} (i.e \textbf{view angle}, \textbf{axial -position}, \textbf{tangential position}) (CTI terminology: \textit{View} -mode) (see Figure \ref{fig:viewgramstorage}) +position}, \textbf{tangential position}) +(see Figure \ref{fig:viewgramstorage}) \end{itemize} This notation means that for a sinogram, the tangential position runs fastest.\\ In both modes, the 3D dataset has been stored in several segments where the number of axial positions in each \textbf{segment} depends -on the \textbf{axial compression (span)}.\\ -Note that we find the CTI terminology confusing, so you will -see it nowhere in STIR except in this document. +on the \textbf{axial compression (span)}. From both modes, 2 different types of 2D datasets can be obtained @@ -175,6 +170,18 @@ \subsection{ view. \end{itemize} +Note that in STIR version 6.0, Time-of-Flight (TOF) will be supported. This introduces another +index. However, \texttt{Sinogram} and \texttt{Viewgram} will remain 2D objects, and \texttt{Segment*} 3D. +This will also be the case once we have layers and energy windows. +In STIR 5.2, we have therefore introduced new classes +\texttt{SinogramIndices}, \texttt{ViewgramIndices}\footnote{Replacing \texttt{ViewSegmentNumbers} in previous versions of STIR.} +and \texttt{SegmentIndices}, containing all ``other'' +indices necessary to get at the corresponding data. This means that the recommended coding style is now +\begin{verbatim} +ViewgramIndices vg_idx(view_num,segment_num); +auto viewgram = proj_data.get_viewgram(vg_idx); +\end{verbatim} +In the next version of STIR, we will introduce extra methods to be able to conveniently loop over all viewgrams etc. \begin{figure}[htbp] \begin{center} @@ -231,8 +238,8 @@ \subsection{ The file \textit{stir}/\textit{common.h} contains general configuration info and tries to iron out some incompatibilities for certain compilers. If you include any STIR \textit{.h} file, you are guaranteed -to have included \textit{stir}/\textit{common.h} as well, hence it should -never be included explicitly by a 'user'. +to have included \textit{stir}/\textit{common.h} as well, hence there is usually no need +to include it explicitly. \subsection{ Namespace} @@ -241,12 +248,9 @@ \subsection{ are within this namespace. The effect of this is that conflicts with other symbols are impossible (except when somebody else is using the same namespace). STIR also uses sub-namespaces for -certain things. This usage will probably be expanded in the future.\\ -For older compilers, the namespace can be disabled (see \textit{stir}/\textit{common.h}). -For this reason, we have introduced some macros such as \textit{START\_NAMESPACE\_STIR.} You -could use the macros in your own code as well, although people -with a compiler that does not support namespace really should -upgrade. +certain things. This usage will probably be expanded in the future. +In most of the code, we use macros such as \textit{START\_NAMESPACE\_STIR} +(originally used for older compilers). \subsection{ Naming conventions } @@ -262,6 +266,7 @@ \subsection{ with \textit{num}, e.g. \textit{num\_gates}. \item the number of an item in a sequence end with \textit{num}, e.g. \textit{gate\_num}. +\footnote{We are slowly starting to use \texttt{\_idx} in a few places as ``indices'' are a more general concept, and avoid confusion between the pre- and postfix use of \texttt{num}.} \item a relative time (normally with respect to the scan start) end with \textit{rel\_time}, e.g. \textit{tracer\_injection\_rel\_time}. @@ -389,15 +394,14 @@ \subsection{ \subsection{\label{ssect:AdvancedCppFeatures} Advanced (?) C++ features used} -We attempted to follow the ANSI C++ standard as close as possible. -We expect to have some marginal problems with a 'strict' ANSI -C++ compiler, although there should now be very few cases (as -gcc warns about a lot of stuff). We have used some preprocessor -macros (see \textit{stir/common.h}) to isolate work-arounds for older +STIR uses C++-11 an dis compatible with C++-20 to the best of our knowledge. +For legacy reasons, we have some preprocessor +macros (see \textit{stir/common.h}) that were used isolate work-arounds for older compilers. Ideally, all these \#ifdefs should disappear at a -later development stage of the library. We gradually give up -on supporting older compilers anyway. +later development stage of the library. +This section describes some C++ features used in STIR that might not be so familiar to +non-C++ programmers. \subsubsection{Templates} The library uses templates very often. This allows us to write @@ -417,23 +421,15 @@ \subsubsection{Templates} actual definition is in a \textit{.txx} file that you can include.\\ We do use partial class template specialisation in some places. However, (very ugly) work-arounds are provided for compilers -that do not support this feature (although not anymore in recent code).\\ -Very occasionally we use member templates (in such a form that -it could be compiled by Visual C++ 6.0, i.e. inline in the class definition). +that do not support this feature (although these areee being gradually removed).\\ +Very occasionally we use member templates. \subsubsection{Run Time Type Information} Another C++ feature that we use is Run Time Type Information (RTTI). We almost exclusively use this to check validity of pointer (or reference) casts down a hierarchy ('down-casting'). See section \ref{sect:classhierarchies} -for an example. If a compiler does not support -RTTI (but all compilers do these days), it would be possible to declare -a (essentially empty) \texttt{dynamic\_cast} -template function such that the above code would compile. The -resulting programme would become inherently type-unsafe though, -and we recommend that you upgrade your compiler. Some code does -rely on explicit type checking, you would have to check this -in detail if you don't have RTTI. +for an example. \subsubsection{ Iterators} @@ -447,7 +443,7 @@ \subsubsection{ \begin{verbatim} void f(vector& v) { - for (vector::iterator iter = v.begin(); iter != v.end(); ++iter) + for (auto iter = v.begin(); iter != v.end(); ++iter) *iter += 2; } \end{verbatim} @@ -464,7 +460,7 @@ \subsubsection{ template void f(Container& v) { - for (typename Container::iterator iter = v.begin(); iter != v.end(); ++iter) + for (auto iter = v.begin(); iter != v.end(); ++iter) *iter += 2; } \end{verbatim} @@ -483,7 +479,7 @@ \subsubsection{ template void f(Array& a) { - for (typename Array::full_iterator iter = a.begin_all(); + for (auto iter = a.begin_all(); iter != a.end_all(); ++iter) *iter += 2; } @@ -511,18 +507,14 @@ \subsubsection{ (something which you cannot achieve with an ordinary pointer). Generally, the destructor of the smart pointer will make sure that the object it points to is deleted (when appropriate).\\ -\textit{std::auto\_ptr} is a standard smart pointer class which is +\textit{std::unique\_ptr} is a standard smart pointer class which is suitable for pointers to objects where there's only one smart -pointer for each object. Unfortunately, its syntax has been decided -rather late in ANSI C++ and its implementation requires some -advanced C++ features, leading to a non-standard implementation -of \texttt{std::auto\_ptr} in older compilers (including VC 6.0) and it is now -superseded by \texttt{std::unique\_ptr} in C++11. For this -reason, we do not use this smart pointer class too much, and +pointer for each object. Unfortunately, a lot of STIR was written +pre-C++-11 and we do not use this smart pointer class too much, and generally use \texttt{shared\_ptr} instead. This is somewhat unfortunate as these two smart pointers are generally quite different concepts. -\texttt{shared\_ptr} is a STIR (or boost) smart pointer class which +\texttt{shared\_ptr} is a wrapper around the \texttt{std} (or boost\footnote{Boost smart pointer are probably no longer supported, and will no longer be from STIR 6.0.}) smart pointer class which is suitable when there are (potentially) more than one pointer pointing to the same object. It keeps a reference count such that the object is (only) deleted when the last shared pointer @@ -561,7 +553,8 @@ \subsubsection{ to test if a pointer is `null', i.e. does not point to anything. \textit{is\_null\_ptr} works with ordinary pointers, shared\_ptrs and auto\_ptrs. Use it to avoid that your code depends on what type of (smart) pointer -you are using. +you are using. Note however that since C++-11, it is more convenient to use the +automatic converstion to \texttt{bool}. \subsection{Generic functionality of STIR classes} This is a (very incomplete) section describing some functionality that @@ -749,8 +742,9 @@ \subsection{ \end{figure} -Currently missing is support for fan-beam data.\\ -List-mode data is experimentally supported since version 1.2. +Currently missing is support for fan-beam data. + +List-mode data is supported as well. See the \textit{Listmode} base-class and its hierarchy. \subsection{ Data (or image) processor hierarchy} @@ -1046,6 +1040,11 @@ \section{ number = 1; a = 2.4F; } + bool post_processing() // will be renamed to post_parsing() + { + // do some checks and handle some extra variables + return false; // everything was ok + } }; int main(int argc, char **argv) @@ -1088,8 +1087,7 @@ \subsection{Images} You will have to include \texttt{stir/IO/write\_to\_file.h} for this to work. This will write using the default output file format and is equivalent to the following: \begin{verbatim} - typedef DiscretisedDensity<3,float> DataType ; - shared_ptr output_format_sptr = + auto output_format_sptr = OutputFileFormat::default_sptr(); output_format_sptr->write_to_file(filename, density); \end{verbatim} @@ -1189,10 +1187,10 @@ \section{ template void f(DiscretisedDensity<3, elemT>& density) { - VoxelsOnCartesianGrid * voxels_ptr = + auto voxels_ptr = dynamic_cast< VoxelsOnCartesianGrid * > (&density); - if (voxels_ptr == NULL) - error("f: can only handle images of type VoxelsOnCartesianGrid\n"); + if (!voxels_ptr) + error("f: can only handle images of type VoxelsOnCartesianGrid"); ... } \end{verbatim} @@ -1286,7 +1284,7 @@ \section{ ) # declare dependencies on other STIR libraries, for instance -target_link_libraries($(dir) buildblock) +target_link_libraries($(dir) PUBLIC buildblock) # add to list of libraries for STIR to include in linking list(APPEND STIR_LIBRARIES $(dir)) diff --git a/documentation/STIR-general-overview.tex b/documentation/STIR-general-overview.tex index ca233c4af..28d8b8626 100644 --- a/documentation/STIR-general-overview.tex +++ b/documentation/STIR-general-overview.tex @@ -17,7 +17,7 @@ \textbf{{\huge STIR \\ General Overview}}\\ \textbf{Kris Thielemans}\\ -\textbf{\textit{version 2.3}} +\textbf{\textit{version 5.2}} \end{spacing} @@ -57,18 +57,21 @@ \section{ Disclaimer} Many names used below are trademarks owned by various companies. They are -is fully acknowledged, but not explicitly stated. +fully acknowledged, but not explicitly stated. + +This document has been brought somewhat up-to-date for version 5.2, but there is more work to do. \section{ Overview} STIR (\textit{Software for Tomographic Image Reconstruction}) is Open Source software (written in C++) consisting of classes, functions -and utilities for 3D PET image reconstruction, although it is +and utilities for 3D PET and SPECT image reconstruction, although it is general enough to accommodate other imaging modalities. An overview of STIR 2.x is given in [Thi12], which you ideally refer to in your paper. +See the STIR website for more details on how to reference STIR, depending on which functionality you use. -STIR consists of 2 parts. +STIR consists of 3 parts. \begin{itemize} \item A library providing building blocks for image and projection @@ -76,6 +79,7 @@ \section{ \item Applications using this library including basic image manipulations, file format conversions and of course image reconstructions. +\item Python interface to the library via SWIG. \end{itemize} The library has been designed so that it can be used for many @@ -85,8 +89,8 @@ \section{ yet. This will enable the software to be run not only on single processors, but also on massively parallel computers, or on clusters of workstations. \\ -STIR is portable on all systems supporting the GNU C++ compiler -or MS Visual C++ (or hopefully any ANSI C++ compliant compiler). +STIR is portable on all systems supporting the GNU C++ compiler, CLang++, Intel C++, +or MS Visual C++ (or hopefully any C++-11 compliant compiler). The library is fully documented.\\ The object-oriented features make this library very modular and flexible. This means that it is relatively easy to add new algorithms, @@ -94,17 +98,16 @@ \section{ It is even possible to select at run-time which version of these components you want to use.\\ The software is \textbf{freely available} for downloading under the -GNU LPGL (the library) or GNU GPL (the applications) license, -see the 'Registration' section of our web-site. \\ +Apache 2.0 license.\\ \textbf{It is the hope of the collaborators of the STIR project that -other researchers in the PET community will use this library +other researchers in the PET and SPECT will use this library for their own work, extending it and making their work available as well.} Please subscribe to some of our mailing lists if you are interested. \\ In its current status, the software is mainly a research tool. It is probably not friendly enough to use in a clinical setting. In any case, \textbf{STIR should not be used to generate images for diagnostic -purposes}, as there is no warranty, and most definitely no FDA approval +purposes}, as there is no warranty, and most definitely no FDA approval nor CE marking @@ -132,10 +135,7 @@ \section{ algorithm type, etc.); \item multi-dimensional arrays (any dimension) with various operations, including numeric manipulations; -\item reading and writing (I/O) data in Interfile format (for -which a 3D PET extension is proposed), reading of GE Advance -sinogram data, limited reading and writing of ECAT6 and ECAT7 -and conversion between ECAT6 and Interfile; +\item reading of various raw data as well as writing in Interfile format; \item classes of projection data (complete data set, segments, sinograms, viewgrams) and images (2D and 3D); \item various filter transfer functions (1D, 2D and 3D); @@ -156,7 +156,7 @@ \section{ \begin{figure}[htbp] \begin{center} \includegraphics[width=5.667in, height=1.137in]{graphics/STIR-general-overviewFig1} -\caption{Current hierarchy for back projectors.} +\caption{Somewhat outdated hierarchy for back projectors.} \end{center} \end{figure} @@ -193,7 +193,7 @@ \subsection{FBP} \subsection{3DRP} The 3DRP [Kin89] algorithm is often considered the 'reference' algorithm for 3D PET. It is a 3D FBP algorithm which uses reprojection -to fill in the missing data. +to fill in the missing data. However, this is not actively maintained anymore. \subsection{ Ordered Subsets Maximum A Posteriori using the One Step Late @@ -203,11 +203,7 @@ \subsection{ its accelerated variant OSEM (Ordered Set Expectation Maximization) [Hud94] are iterative methods for computing the maximum likelihood estimate of the tracer distribution based on the measured projection -data. Their benefits have been examined in various reports. Notably, -it has been found (e.g. [Rea98b]) that 3D OSEM can produce images -of higher quality, in terms of both resolution and contrast, -than analytic algorithms of common commercial use such as the -reprojection (3DRP) algorithm [Kin89]. +data. One drawback of OSEM is its tendency to develop noise artefacts with increasing iterations. As a remedy, various modifications @@ -371,47 +367,17 @@ \section{ \section{ Currently supported systems} - -These are the systems we tested our software on. We think that -it will work in other cases as well with minor tuning. Experience -will tell\dots - - - -\subsection{ -Computing platforms} +We regularly run STIR on Lniux, MacOS and Windows. Check our GitHub Actions +and Appveyor for details. -\subsubsection{ -Serial versions of algorithms} -A lot of STIR was developed kust before 2000. Therefore, the code -contains in various places work-arounds for deficient C++ compilers -(e.g. GNU g++ 2.95.2 or Microsoft Visual C++ 6.0 service pack 3). -However, we are no longer -testing STIR on such older systems. It might be possible to revive this, -but we strongly recommend that you get a recent compiler instead. -\begin{itemize} -\item -Intel based PCs using Windows XP or higher (NT probably still works if you -find a compiler for it). These were tested with Microsoft Visual C++ 2003, 2005 up to -Visual Studio 2010, -and various GNU g++ versions (up to at least 4.5.3) via CYGWIN (http://cygwin.com). -\item -Linux using GNU g++ (from 3.5 to 4.6). -Intel's C++ compiler 12.1.0 works fine (tested on Linux), older versions had some problems with the Fourier transform code but the reconstruction itself was fine. -\item -MacOS using GNU g++. -\end{itemize} -However, people run STIR on various other systems, including Solaris. Please -check the mailing lists. - \textbf{Warning}: We currently have a problem in the incremental backprojection routines due to different rounding -of floating point calculations on Sparc, HP and Opteron processors. -Thanks to a fix of Bing Bai this problem might have disappeared on your processor, but it is still there on recent 64-bit processors. +of floating point calculations. You will find out if this problem still exists when you run the recon\_text\_pack available on the STIR web-site. Please let us know. +This currently only affects the FBP routines. See the User's Guide for how to use another backprojector. @@ -419,7 +385,8 @@ \subsubsection{ Parallel versions of algorithms} \begin{itemize} \item a parallel version of \texttt{OSMAPOSL} and \texttt{OSSPS} using MPI -\item a threaded version of FBP2D using OPENMP, but for most systems, there is no performance gain. +\item OpenMP versions of many functions +\item CUDA versions for NiftyPET and parallelproj \end{itemize} @@ -427,18 +394,7 @@ \subsubsection{ \subsection{ PET scanners} -\begin{itemize} -\item -GE Advance, Discovery LS, Discovery ST, Discvoery STE, Discovery RX, -Discovery 600 -\item -PRT-1 -\item -ECAT 921, 931, 953, 951, 962, 966 (and most other ECAT systems) -\item -HRRT apparently works as well -\item Philips Allegro (using its convention of having a ``virtual'' extra crystal to cover the gaps). -\end{itemize} +See \texttt{buildblock/Scanners.cxx}. Other cylindrical PET scanners can easily be added. This can be done without modifying the code (see the Wiki). @@ -449,37 +405,27 @@ \subsection{ \begin{itemize} \item An extension of the Interfile standard (see the ''Other info'' section of the STIR website.) +\item +ITK can be used to read many image file formats +\item Siemens ``interfile-like'' data +\item GE RDF9 (if you have the HDF5 libraries) \item The commercial AnalyzeAVW library from the BIR, Mayo University, can be used -to read images via a conversion utility. +to read images via a conversion utility. This will be dropped in version 6.0. \item -GE Advance VOLPET sinogram format, but for reading only (with -very limited support for the header). If you have a research contract with GE, -you might be able to get a STIR extension to read Discovery ST, STE and RX data in both -VOLPET and RDF format. +GE Advance VOLPET sinogram format, but for reading only. This will be dropped in version 6.0. \item -ECAT 7 matrix format for reading only (for reading multiple -frames/gates, an Interfile -header has to be created using a utility we provide). Single -frame/gate/bed images can be directly read or written. Single -frame/gate/bed sinograms can be directly read. +ECAT 7 matrix format for reading only might work but is no longer supported. However, this file format needs the ECAT Matrix library (developed previously by M. Sibomana and C. Michel at Louvain la Neuve, Belgium). This library is no longer maintained however. -\item -ECAT 6 matrix format was supported in previous version of STIR, but -is currently broken on some systems (probably due to a bug in the ECAT -MATRIX library for ECAT 6). ECAT6 was supported by using conversion -utilities to and from -Interfile. Single frame/gate/bed images could be directly read -or written. \end{itemize} See the User's Guide for more detail on the supported file formats.\\ In addition, a separate set of classes is avalailable to read list-mode data. Only a few scanners are currently supported (such as the -ECAT HR+ and HR++), although it should not be too difficult to +ECAT HR+ and HR++, GE RDF9, Siemens mMR and SAFIR), although it should not be too difficult to add your own (if you know the list-mode file format!). diff --git a/documentation/history.htm b/documentation/history.htm index 076f96d3b..c4a258e8b 100644 --- a/documentation/history.htm +++ b/documentation/history.htm @@ -11,6 +11,8 @@

History of public releases of the STIR software

The following links give you a summary what has changed. However, the ChangeLog is the definite (but tedious) information.
    +
  • version 5.2.0 (dated 30 October 2023)
    • +
  • version 5.1.2 (dated 10 September 2023)
  • version 5.1.1 (dated 26 Augustus 2023)
  • version 5.1.0 (dated 14 Jan 2023)
  • version 5.0.2 (dated 26 May 2022)
    • @@ -46,7 +48,7 @@

    History of public releases of the PARAPET software

    -Last modified: January 02 2023 +Last modified: 29 October 2023 diff --git a/documentation/release_5.2.htm b/documentation/release_5.2.htm index 53462b803..315acfee8 100644 --- a/documentation/release_5.2.htm +++ b/documentation/release_5.2.htm @@ -7,7 +7,7 @@

    Summary of changes in STIR release 5.2

    -

    This version is 100% backwards compatible with STIR 5.0. However, there is a change in the output of scatter estimation, see below for more information. +

    This version is 100% backwards compatible with STIR 5.0 as far as usage goes. However, there are changes in the output of scatter estimation and ECAT8 normalisation, see below for more information.

    Overall summary

    @@ -15,11 +15,12 @@

    Overall summary

    improvements to the documentation. See also the 5.2 milestone on GitHub.

    -

    Overall code management and assistance by Kris Thielemans (UCL and ASC).

    +

    Overall code management and assistance by Kris Thielemans (UCL and ASC). Other main contributors + were Daniel Deidda (NPL) and Markus Jehl (Positrigo).

    Patch release info

      -
    • 5.2.0 released ??/??/2023
      • +
    • 5.2.0 released 30/10/2023

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

    @@ -30,6 +31,27 @@

    Bug fixes

  • Setting SPECTUB resolution model with STIR python or SIRF divided slope by 10 in error. The problem did not occur when set using parameter file
+

Changed functionality

+
    +
  • The ECAT8 normalisation (used for the Siemens mMR) code now takes the 4th component axial effects into account. + These normalisation factors are therefore different (even up to ~10%). This gives improved axial uniformity in the images. + The use of the axial effects can be switched off by adding setting use_axial_effects_factors:=0 to the + parameter file (see an example in examples/Siemens-mMR/correct_projdata_no_axial_effects.par), or the class member + of the same name.
    + In addition, the Siemens normalisation header is now read (using a new class InterfileNormHeaderSiemens) such that + hard-coded variables for the Siemens mMR have been removed. Further testing of this functionality is still required however. +
    PR #1182. +
    • +
  • Interfile header parsing now correctly identifies keywords that contain a colon by checking for :=.
    • +
  • + The set_up() method of the ray-tracing projection matrix now skips further processing + if it was already called with data of the same characteristics. This will means that any cached data + will be re-used, potentially leading to a speed-up when re-using it from Python. +
    PR #1281. +
    • + +
+

New functionality

  • The Discretised Shape3D shape/ROI has now an extra value label index. For ROIs, this allows @@ -61,16 +83,33 @@

    New functionality

    but you could use it to optimise the number of OpenMP threads to use for your data.
    PR #1237.
    • +
  • New classes SegmentIndices, ViewgramIndices and SinogramIndices, used by ProjData related classes, as opposed + to having to specify all the elements directly, e.g. in C++ + 
    +      auto sinogram = proj_data.get_sinogram(sinogram_indices);
+
    + This makes these functions more future proof, in particular for TOF. The older functions are now deprecated. + Note that as Bin is now derived from ViewgramIndices, instantations of Bin can now be used to specify the indices as well in most places. +
    + There is still more work to do here, mostly related to the symmetries. +
    PR #1273. +
-

New examples

+

Python (and MATLAB)

    +
  • Examples use stir.ProjData.read_from_file as opposed to stir.ProjData_read_from_file. The former is supported since SWIG 3.0, and the default from SWIG 4.1. +
    • +
  • Addition of DetectionPosition and DetectionPositionPair.
    • +
  • bin.time_frame_num is now no longer a function in Python, but acts like a variable + (as the other Bin members). +
    • +
  • Addition of RadionuclideDB
-

Python

+

New examples

    -
  • Examples use stir.ProjData.read_from_file as opposed to stir.ProjData_read_from_file. The former is supported since SWIG 3.0, and the default from SWIG 4.1. -
    • +
  • examples/python/construct_projdata_demo.py illustrates constructing a ProjDataInMemory

Changed functionality

@@ -95,18 +134,30 @@

Changed functionality

This is to reduce the default output during iterative reconstructions.
PR #1243. +
  • The Succeeded class has a new method bool succeeded() enabling more concise code (avoiding the need for comparing with Succeeded::yes which is especially verbose in Python). +
    • +
  • The example files for the Siemens mMR now use lower min/max thresholds for the (single) scatter scale. This gives better results, see Issue #1163. +
    PR #1279. +
    • -

    Deprecated functionality

    +

    Deprecated functionality and upcoming changes to required tool versions

    • The following functions (previously used for upsampling the scatter estimate) have been made obsolete or replaced, and will be removed in STIR version 6.0.0: interpolate_axial_position, extend_sinogram_in_views and extend_segment_in_views
      • +
    • Constructors/functions in ProjData related classes that explicitly use axial_pos_num, view_num etc in their arguments are now deprecated, + and should be replaced by their respective versions that use SegmentIndices, ViewgramIndices or SinogramIndices. The former will not be + compatible with TOF information that will be introduced in version 6.0.0. +
      • +
    • Use of the AVW library to read Analyze files will be removed in 6.0, as this has not been checked in more than 15 years. Use ITK instead.
      • +
    • GE VOLPET and IE support will be removed in 6.0, as we have no files to test this, and it's obsolete anyway.
      • +
    • STIR version 6.0.0 will require C++ 14 (currently we require C++ 11, but already support C++ 20) and CMake 3.14.

    Build system and dependencies

      -
    • CMake 3.12 is now required on Windows. (It will be required on all platforms in the next version).
      • +
    • CMake 3.12 is now required on Windows.
    • We now use CMake's OBJECT library feature for the registries. This avoids re-compilation of the registries for every executable and therefore speeds-up building time. Use of STIR in an external project is not affected as long as the recommended practice was followed. This is now documented in the User's Guide.
      PR #1141.
      • @@ -158,10 +209,12 @@

      Minor (?) bug fixes

      Documentation changes

        +
      • Updated the STIR developers guide, which was quite out-of-date w.r.t. C++ features etc.

      recon_test_pack changes

        +
      • Updated headers of most images and projection data to avoid warnings.

      Other changes to tests

      diff --git a/documentation/release_notes_TOF.htm b/documentation/release_notes_TOF.htm index 96aadb367..068ccd28f 100644 --- a/documentation/release_notes_TOF.htm +++ b/documentation/release_notes_TOF.htm @@ -7,12 +7,12 @@

      Summary of changes in STIR release 6.0

      -

      This version is 95% backwards compatible with STIR 4.0 for the user (see below). +

      This version is 95% backwards compatible with STIR 5.x for the user (see below). Developers might need to make code changes as detailed below.

      Overall summary

      -

      +

      This release is a major upgrade adding Time of Flight (TOF) capabilities to STIR.

      @@ -58,6 +58,8 @@

      New functionality

      Changed functionality

      • + We now always check (ProjDataInfo*NoArcCorr) if number of tangential positions in the projection data exceeds the maximum number + of non arc-corrected bins set for the scanner. If it is, an error is raised.
      @@ -77,7 +79,7 @@

      Known problems

      Minor bug fixes

        -
      • +
      diff --git a/examples/PET_simulation/scatter_template.hs b/examples/PET_simulation/scatter_template.hs index 9929fc413..e8086c80f 100644 --- a/examples/PET_simulation/scatter_template.hs +++ b/examples/PET_simulation/scatter_template.hs @@ -32,10 +32,10 @@ Number of detectors per ring := 64 Inner ring diameter (cm) := 88.62 Average depth of interaction (cm) := 0.84 Distance between rings (cm) := 1.962 -Default bin size (cm) := -1 +Default bin size (cm) := -1 ; unused as not arc-correcting View offset (degrees) := 0 -Maximum number of non-arc-corrected bins := 0 -Default number of arc-corrected bins := 0 +Maximum number of non-arc-corrected bins := 35 +Default number of arc-corrected bins := 0 ; unused as not arc-correcting Number of blocks per bucket in transaxial direction := 0 Number of blocks per bucket in axial direction := 0 Number of crystals per block in axial direction := 0 diff --git a/examples/ROOT_files/ROOT_STIR_consistency/.gitignore b/examples/ROOT_files/ROOT_STIR_consistency/.gitignore new file mode 100644 index 000000000..72b44edd9 --- /dev/null +++ b/examples/ROOT_files/ROOT_STIR_consistency/.gitignore @@ -0,0 +1 @@ +pretest_output/ diff --git a/examples/ROOT_files/ROOT_STIR_consistency/DebugScripts/debug_consistency_with_root.py b/examples/ROOT_files/ROOT_STIR_consistency/DebugScripts/debug_consistency_with_root.py index 82b0cb163..ef06e7395 100644 --- a/examples/ROOT_files/ROOT_STIR_consistency/DebugScripts/debug_consistency_with_root.py +++ b/examples/ROOT_files/ROOT_STIR_consistency/DebugScripts/debug_consistency_with_root.py @@ -87,15 +87,15 @@ def point_cloud_3D(data_handler): ax = fig.add_subplot(projection='3d') # Plot all the points (intensity increases for multiple points) - ax.scatter(data_handler.voxel_coords[:, 0], data_handler.voxel_coords[:, 1], data_handler.voxel_coords[:, 2]) + ax.scatter(data_handler.voxel_coords[:, 0], data_handler.voxel_coords[:, 1], data_handler.voxel_coords[:, 2], label='Most Likely Positions') # Plot the original point and tolerance ox = data_handler.original_coord[0] oy = data_handler.original_coord[1] oz = data_handler.original_coord[2] tol = data_handler.tolerance - ax.plot(ox, oy, oz, c='r', marker='o') - # plot tolerence around original point + ax.plot(ox, oy, oz, c='r', marker='o', label='Origin and Tolerance') + # plot tolerance around original point ax.plot([ox + tol, ox - tol], [oy, oy], [oz, oz], c='r', marker="_", label='_nolegend_') ax.plot([ox, ox], [oy + tol, oy - tol], [oz, oz], c='r', marker="_", label='_nolegend_') ax.plot([ox, ox], [oy, oy], [oz + tol, oz - tol], c='r', marker="_", label='_nolegend_') @@ -107,7 +107,7 @@ def point_cloud_3D(data_handler): xerror = np.std(data_handler.voxel_coords[:, 0]) yerror = np.std(data_handler.voxel_coords[:, 1]) zerror = np.std(data_handler.voxel_coords[:, 2]) - ax.plot(fx, fy, fz, linestyle="None", marker="o", c='g') + ax.plot(fx, fy, fz, linestyle="None", marker="o", c='g', label='Mean coords and stddev') ax.plot([fx + xerror, fx - xerror], [fy, fy], [fz, fz], marker="_", c='g', label='_nolegend_') ax.plot([fx, fx], [fy + yerror, fy - yerror], [fz, fz], marker="_", c='g', label='_nolegend_') ax.plot([fx, fx], [fy, fy], [fz + zerror, fz - zerror], marker="_", c='g', label='_nolegend_') @@ -115,7 +115,7 @@ def point_cloud_3D(data_handler): ax.set_xlabel('x (mm)') ax.set_ylabel('y (mm)') ax.set_zlabel('z (mm)') - ax.legend(['Voxel Positions', 'Origin and Tolerance', 'Mean coords and stddev']) + ax.legend() plt.show() @@ -173,7 +173,7 @@ def point_cloud_3D_all(dict_of_data_handlers): plt.show() -def __extract_data_from_csv_file(filename): +def extract_data_from_csv_file(filename): """ Loads data from a csv file and returns a numpy array containing the data. Assumes first column is the row's key. @@ -222,7 +222,7 @@ def __init__(self, filename): print(f"Loading data: {self.filename}") # Extract the original coordinate and voxel coordinates as 2D numpy arrays - self.original_coord, self.voxel_coords, self.tolerance = __extract_data_from_csv_file(filename) + self.original_coord, self.voxel_coords, self.tolerance = extract_data_from_csv_file(filename) if self.tolerance is None: raise Exception("No tolerance information found in file") @@ -387,8 +387,8 @@ def TOF_evaluation(filename_prefix, file_extension=".csv"): # Main Script # ===================================================================================================== def main(): - print("\nUSAGE: After `make test` or `test_view_offset_GATE` has been run,\n" - "run `debug_view_offset_consistency` from `ROOT_STIR_consistency` directory or input that directory as an " + print("\nUSAGE: After `make test` or `test_consistency_with_GATE` has been run,\n" + "run `debug_consistency_with_root.py` from `ROOT_STIR_consistency` directory or input that directory as an " "argument.\n") # Optional argument to set the directory of the output of the test_consistency_with_root CTest. diff --git a/examples/ROOT_files/ROOT_STIR_consistency/README.md b/examples/ROOT_files/ROOT_STIR_consistency/README.md index 476fa0f7d..e2798f368 100644 --- a/examples/ROOT_files/ROOT_STIR_consistency/README.md +++ b/examples/ROOT_files/ROOT_STIR_consistency/README.md @@ -8,14 +8,17 @@ See STIR/LICENSE.txt for details These files were included to : * Test non-TOF ROOT and STIR consistency, particularly the rotation. -* Test the TOF STIR implementation was correct (NOT YET IMPLEMENTED). +* Test the TOF STIR implementation is correct. + +See src/recon_test/test_consistency_with_GATE.cxx. This test is run +automatically when using ctest. Directories ----- -- `SourceFiles/`: Contains 8 `generate_image` parameter files and GATE macro files for the emission source positions. One pair of files for each test. +- `SourceFiles/`: Contains 8 `generate_image` parameter files and GATE macro files for the emission source positions. One pair of files for each simulation. - `Gate_macros/`: Contains the GATE macro files for generating the data. -- `DebugScripts`: Contains scripts for better understanding the tests. +- `DebugScripts/`: Contains scripts for better understanding the tests. FILES @@ -30,14 +33,15 @@ ______ Methodology ---- - 1. Generate the ROOT data. - 1. Run `./run_pretest_script.sh` in the terminal to generate the ROOT files (requires Gate) for different point sources, or - 2. Download the ROOT data and proceed without Gate simulation. + 1. Get the ROOT data: either + - Run `./run_pretest_script.sh` in the terminal to generate the ROOT files (requires Gate) for different point sources, or + - Download the ROOT data and proceed without Gate simulation. This is done + by ctest when the STIR build was configured with `DOWNLOAD_ZENODO_TEST_DATA=ON`. - 2. Run the STIR test: `src/recon_test/test_view_offset_root`. + 2. Run the STIR test: `src/recon_test/test_consistency_with_GATE` (Best via ctest). This test should tell you whether it failed or not by testing if the LOR passes by, or close to, the original point source position. - 3. Run the python scripts in `DebugScripts` to better understand erros and to give a more in-depth analysis. + 3. Run the python scripts in `DebugScripts` to better understand errors and to give a more in-depth analysis. _____ diff --git a/examples/Siemens-mMR/README.md b/examples/Siemens-mMR/README.md index 85dc711cd..40621c019 100644 --- a/examples/Siemens-mMR/README.md +++ b/examples/Siemens-mMR/README.md @@ -16,9 +16,9 @@ The scripts don't do proper error handling. If they suddenly stop without diagno and look for a log file. The scripts need template files from this directory. By default -they assume they are located in `~/devel/STIR/examples/Siemens-mMR`. -If that is not the case, you can set the pardir variable as above, or -by doing for instance +they assume they are located in the `Siemens-mMR` directory of the STIR +examples folder, as obtained by `stir_config --examples-dir`. +This can be overridden by doing for instance ```sh pardir=~/STIR/examples/Siemens-mMR diff --git a/examples/Siemens-mMR/correct_projdata.par b/examples/Siemens-mMR/correct_projdata.par index a678cce16..246edccad 100644 --- a/examples/Siemens-mMR/correct_projdata.par +++ b/examples/Siemens-mMR/correct_projdata.par @@ -40,6 +40,9 @@ correct_projdata Parameters := Bin Normalisation type := from ecat8 Bin Normalisation From ecat8 := normalisation filename:= ${ECATNORM} + ; keyword that can be used to write the components to a separate text files for debugging + ; files are written in the current directory and are called geom_out.txt etc. + ; write_components_to_file := 0 End Bin Normalisation From ecat8:= ; attenuation image, will be forward projected to get attenuation factors diff --git a/examples/Siemens-mMR/correct_projdata_no_axial_effects.par b/examples/Siemens-mMR/correct_projdata_no_axial_effects.par new file mode 100644 index 000000000..861341516 --- /dev/null +++ b/examples/Siemens-mMR/correct_projdata_no_axial_effects.par @@ -0,0 +1,62 @@ +correct_projdata Parameters := + ; a sample file that switches the use of the "axial effects" + ; (the 4th component in the ECAT8 norm file) off + input file := ${INPUT} + + ; Current way of specifying time frames, pending modifications to + ; STIR to read time info from the headers + ; see class documentation for TimeFrameDefinitions for the format of this file + ; time frame definition filename := frames.fdef + + ; if a frame definition file is specified, you can say that the input data + ; corresponds to a specific time frame + ; the number should be between 1 and num_frames and defaults to 1 + ; this is currently only used to pass the relevant time to the normalisation + ; time frame number := 1 + + ; output file + ; for future compatibility, do not use an extension in the name of the + ; output file. It will be added automatically + output filename := ${OUTPUT} + + ; default value for next is -1, meaning 'all segments' + ; maximum absolute segment number to process := + + + ; use data in the input file, or substitute data with all 1's + ; (useful to get correction factors only) + ; default is '1' + use data (1) or set to one (0) := 0 + + ; precorrect data, or undo precorrection + ; default is '1' + ; apply (1) or undo (0) correction := + + ; parameters specifying correction factors + ; if no value is given, the corresponding correction will not be performed + + ; random coincidences estimate, subtracted before anything else is done + ;randoms projdata filename := random.hs + ; normalisation (or binwise multiplication, so can contain attenuation factors as well) + Bin Normalisation type := from ecat8 + Bin Normalisation From ecat8 := + normalisation filename:= ${ECATNORM} + use_axial_effects_factors:=0 + End Bin Normalisation From ecat8:= + + ; attenuation image, will be forward projected to get attenuation factors + ; OBSOLETE + ;attenuation image filename := attenuation_image.hv + + ; forward projector used to estimate attenuation factors, defaults to Ray Tracing + ; OBSOLETE + ;forward_projector type := Ray Tracing + + ; scatter term to be subtracted AFTER norm+atten correction + ; defaults to 0 + ;scatter projdata filename := scatter.hs + + ; to interpolate to uniform sampling in 's', set value to 1 + ; arc correction := 0 + +END:= diff --git a/examples/Siemens-mMR/correct_projdata_only_counts.par b/examples/Siemens-mMR/correct_projdata_only_counts.par index 734ba63b9..68cb0ff0b 100644 --- a/examples/Siemens-mMR/correct_projdata_only_counts.par +++ b/examples/Siemens-mMR/correct_projdata_only_counts.par @@ -1,4 +1,6 @@ correct_projdata Parameters := +; a par file to investigate what BinNormalisationFromECAT8 does when only +; taking span/mashing into account (i.e. actually ignoring the norm factors themselves) input file := ${INPUT} @@ -44,7 +46,8 @@ correct_projdata Parameters := use_detector_efficiencies:=0 use_geometric_factors:=0 use_crystal_interference_factors:=0 - End Bin Normalisation From ecat8:= + use_axial_effects_factors:=0 + End Bin Normalisation From ecat8:= ; attenuation image, will be forward projected to get attenuation factors ; OBSOLETE diff --git a/examples/Siemens-mMR/scatter_and_recon.sh b/examples/Siemens-mMR/scatter_and_recon.sh index 9442840da..fa45238f3 100755 --- a/examples/Siemens-mMR/scatter_and_recon.sh +++ b/examples/Siemens-mMR/scatter_and_recon.sh @@ -75,7 +75,7 @@ echo "Estimating scatter (be patient). Log saved in output/scatter.log" # filename-prefix for additive sino (i.e. "precorrected" sum of scatter and randoms) total_additive_prefix=output/total_additive num_scat_iters=3 -scatter_pardir=${pardir}/../samples/scatter_estimation_par_files +scatter_pardir=${pardir}/scatter_estimation_par_files # you might have to change this for a different scanner than the mMR scatter_recon_num_subiterations=21 scatter_recon_num_subsets=21 diff --git a/examples/Siemens-mMR/scatter_estimation_par_files/README.md b/examples/Siemens-mMR/scatter_estimation_par_files/README.md new file mode 100644 index 000000000..1bd540f3b --- /dev/null +++ b/examples/Siemens-mMR/scatter_estimation_par_files/README.md @@ -0,0 +1,19 @@ +# Example files for running scatter estimation for the Siemens mMR + +Files made by Nikos Efthimou and fine-tuned by Kris Thielemans.
      +Copyright University of Hull 2018-2019
      +copyright University College London 2016, 2020
      +Distributed under the Apache 2.0 License + +These files are almost identical to those in +[examples/samples/scatter_estimation_par_files/](../../samples/scatter_estimation_par_files/README.md), +see there for some more information. + +Currently the only difference are the lower values for +``` +maximum scatter scaling factor := .5 +minimum scatter scaling factor := 0.1 +``` + +These have been shown to work better for mMR data, see e.g. +[STIR issue #1163](https://github.com/UCL/STIR/issues/1163). diff --git a/examples/Siemens-mMR/scatter_estimation_par_files/postfilter_Gaussian_for_mask.par b/examples/Siemens-mMR/scatter_estimation_par_files/postfilter_Gaussian_for_mask.par new file mode 100644 index 000000000..d6af00a35 --- /dev/null +++ b/examples/Siemens-mMR/scatter_estimation_par_files/postfilter_Gaussian_for_mask.par @@ -0,0 +1,12 @@ +PostFilteringParameters := + Postfilter type := Separable Gaussian +Separable Gaussian Filter Parameters := +x-dir filter FWHM (in mm):= 20 +y-dir filter FWHM (in mm):= 20 +z-dir filter FWHM (in mm):= 15 +; optionally restrict kernel sizes +; x-dir maximum kernel size := 129 +; y-dir maximum kernel size := 129 +; z-dir maximum kernel size := 31 +END Separable Gaussian Filter Parameters := +End PostFiltering Parameters:= diff --git a/examples/Siemens-mMR/scatter_estimation_par_files/run_reconstruction.par b/examples/Siemens-mMR/scatter_estimation_par_files/run_reconstruction.par new file mode 100644 index 000000000..166bdff2e --- /dev/null +++ b/examples/Siemens-mMR/scatter_estimation_par_files/run_reconstruction.par @@ -0,0 +1,44 @@ +Reconstruction Parameters := +reconstruction type := OSMAPOSL +OSMAPOSLParameters := + +objective function type:= PoissonLogLikelihoodWithLinearModelForMeanAndProjData +PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:= +maximum absolute segment number to process := -1 + +projector pair type := Matrix + Projector Pair Using Matrix Parameters := + Matrix type := Ray Tracing + Ray tracing matrix parameters := + number of rays in tangential direction to trace for each bin:= 5 + End Ray tracing matrix parameters := + End Projector Pair Using Matrix Parameters := + +;recompute sensitivity := 0 +;subset sensitivity filenames := scatter_subset_sens_%d.hv + +; reconstruct at large voxel size to save time +zoom := 0.2 + +end PoissonLogLikelihoodWithLinearModelForMeanAndProjData Parameters:= + +; initial estimate := +enforce initial positivity condition:=1 + +number of subsets:= ${scatter_recon_num_subsets} +number of subiterations:=${scatter_recon_num_subiterations} +;save estimates at subiteration intervals:= ${scatter_recon_num_subiterations} + +; smooth a bit as we use a down-sampled scanner (during the scatter estimation resolution can be low) +post-filter type := Separable Gaussian +Separable Gaussian Filter Parameters := + x-dir filter FWHM (in mm):= 15 + y-dir filter FWHM (in mm):= 15 + z-dir filter FWHM (in mm):= 15 +END Separable Gaussian Filter Parameters := +; +; Disable output +; +disable output := 1 +End OSMAPOSLParameters:= +End reconstruction Parameters:= diff --git a/examples/Siemens-mMR/scatter_estimation_par_files/scatter_estimation.par b/examples/Siemens-mMR/scatter_estimation_par_files/scatter_estimation.par new file mode 100644 index 000000000..771f0daf8 --- /dev/null +++ b/examples/Siemens-mMR/scatter_estimation_par_files/scatter_estimation.par @@ -0,0 +1,86 @@ +Scatter Estimation Parameters := + +;Run in debug mode +;; A new folder called extras will be created, in which many +;; extra files will be stored +run in debug mode := 1 + +; Measured data +input file := ${sino_input} + +; Attenuation Image +attenuation image filename := ${atnimg} + +; Normalisation coefficients & attenuation data +Normalisation type := from ProjData + Bin Normalisation From ProjData := + normalisation projdata filename:= ${NORM} + End Bin Normalisation From ProjData:= + +attenuation correction factors filename := ${acf3d} + +;; Background data (not normalised). +; Should be set to the randoms estimate (unless you precorrected, but we haven't tested that) +background projdata filename := ${randoms3d} + +; Mask for tail-fitting +; It will be computed by masking the attenuation image, and forward projecting that. +; If !recompute mask projdata then the filename must be set. +recompute mask projdata := 1 +mask projdata filename := ${mask_projdata_filename} + +; Input or output filename - depends on recompute +recompute mask image := 1 +mask image filename := ${mask_image} +; threshold to be applied after filtering (in cm^-1). Default value is below +mask attenuation image min threshold := 0.003 +; optional filename to specify a filter before thresholding the attenuation image +; By default a Gaussian filter with FWHM (15,20,20) will be used. Here we use an explicit file as an example. +mask attenuation image filter filename := ${scatter_pardir}/postfilter_Gaussian_for_mask.par +;End of Mask + +;Parameter file for the tail fitting of the scatter data (within the mask) +tail fitting parameter filename := ${scatter_pardir}/tail_fitting.par + +; Run simulation and reconstruction in 2D and export SSRB sinograms (currently required) +run in 2d projdata := 1 + +; ScatterSimulation parameters +; could read from a file, but instead we have them below +; scatter simulation parameter filename := ${scatter_pardir}/scatter_simulation.par +Scatter Simulation type := PET Single Scatter Simulation + PET Single Scatter Simulation Parameters := + ; could change some parameters here if you need to (not recommended) + End PET Single Scatter Simulation Parameters:= + +; next option is the default +use scanner downsampling in scatter simulation := 1 + +; could add parameters below, but reading it from file +; reconstruction type := ... +reconstruction parameter filename := ${scatter_pardir}/run_reconstruction.par + +; +; This is the number of times which the Scatter Estimation will +; iterate. Default is 5 + +number of scatter iterations := ${num_scat_iters} + +; Average the first two activity images +do average at 2 := 1 + +; Export scatter estimates of each iteration +export scatter estimates of each iteration := 1 + +output scatter estimate name prefix := ${scatter_prefix} +output additive estimate name prefix:= ${total_additive_prefix} + +maximum scatter scaling factor := 0.4 +minimum scatter scaling factor := 0.1 + +;Upsample and fit +; defaults to 3. +upsampling half filter width := 3 +remove interleaving before upsampling := 1 + +End Scatter Estimation Parameters := diff --git a/examples/Siemens-mMR/scatter_estimation_par_files/tail_fitting.par b/examples/Siemens-mMR/scatter_estimation_par_files/tail_fitting.par new file mode 100644 index 000000000..dc0d1459f --- /dev/null +++ b/examples/Siemens-mMR/scatter_estimation_par_files/tail_fitting.par @@ -0,0 +1,4 @@ +CreateTailMaskFromACFs := + ACF-threshold := 1.1 + safety-margin := 4 +END CreateTailMaskFromACFs := diff --git a/examples/python/construct_projdata_demo.py b/examples/python/construct_projdata_demo.py new file mode 100644 index 000000000..791f22256 --- /dev/null +++ b/examples/python/construct_projdata_demo.py @@ -0,0 +1,47 @@ +# Demo of how to use STIR from Python to construct some projection data from scratch + +# Copyright 2023 - University College London +# This file is part of STIR. +# +# SPDX-License-Identifier: Apache-2.0 +# +# See STIR/LICENSE.txt for details + +#%% Initial imports +import stir + +#%% check list of predefined scanners +print(stir.Scanner.get_names_of_predefined_scanners()) + +#%% create a scanner +scanner=stir.Scanner.get_scanner_from_name("Siemens mMR") + +#%% create a ProjDataInfo, describing the geometry of the data acquired for that scanner +span = 11 +max_ring_diff = 60 +num_views = scanner.get_num_detectors_per_ring()//2 +num_tangential_poss = scanner.get_default_num_arccorrected_bins() +proj_data_info = stir.ProjDataInfo.construct_proj_data_info(scanner, span, max_ring_diff, num_views, num_tangential_poss, False); + +#%% create radionuclide +modality = stir.ImagingModality(stir.ImagingModality.PT) +db = stir.RadionuclideDB() +r = db.get_radionuclide(modality, "^18^Fluorine") + +#%% supported values are determined by your radionuclide JSON database +with open(stir.find_STIR_config_file("radionuclide_names.json")) as f: + print(f.read()) + +#%% create ExamInfo object +exam_info = stir.ExamInfo(modality) +exam_info.set_radionuclide(r) +exam_info.patient_position=stir.PatientPosition(stir.PatientPosition.HFS) +print(exam_info.parameter_info()) + +#%% create projection data in memory, filled with 1 +proj_data = stir.ProjDataInMemory(exam_info, proj_data_info) +proj_data.fill(1) + +#%% write it to file +# STIR can currently only write Interfile projection data +proj_data.write_to_file('ones.hs') diff --git a/examples/python/recon_demo.py b/examples/python/recon_demo.py index 1a65ecbd3..815190573 100755 --- a/examples/python/recon_demo.py +++ b/examples/python/recon_demo.py @@ -48,7 +48,7 @@ print("Enabling interactive-mode for plotting failed. Continuing.") s = recon.set_up(target) -if (s == stir.Succeeded(stir.Succeeded.yes)): +if (s.succeeded()): pylab.figure() for iter in range(1,num_subiterations+1): print('\n--------------------- Subiteration ', iter) diff --git a/recon_test_pack/RPTsens_seg3_PM.ahv b/recon_test_pack/RPTsens_seg3_PM.ahv deleted file mode 100644 index 5e56e67da..000000000 --- a/recon_test_pack/RPTsens_seg3_PM.ahv +++ /dev/null @@ -1,18 +0,0 @@ -!INTERFILE := -!name of data file := RPTsens_seg3_PM.v -!total number of images := 31 -!data offset in bytes := 0 -!imagedata byte order := LITTLEENDIAN -!number format := short float -!number of bytes per pixel := 4 -matrix axis label [1] := x -!matrix size [1] := 60 -scaling factor (mm/pixel) [1] := 4.44114 -matrix axis label [2] := y -!matrix size [2] := 60 -scaling factor (mm/pixel) [2] := 4.44114 -;Correct value is of keyword (commented out) -;!slice thickness (pixels) := 0.759939 -;Value for Analyze -!slice thickness (pixels) := 3.375 -!END OF INTERFILE := diff --git a/recon_test_pack/RPTsens_seg4.ahv b/recon_test_pack/RPTsens_seg4.ahv deleted file mode 100644 index 28ff3cb02..000000000 --- a/recon_test_pack/RPTsens_seg4.ahv +++ /dev/null @@ -1,15 +0,0 @@ -!INTERFILE := -!name of data file := RPTsens_seg4.v -!total number of images := 31 -!data offset in bytes := 0 -!imagedata byte order := LITTLEENDIAN -!number format := float -!number of bytes per pixel := 4 -matrix axis label [1] := x -!matrix size [1] := 97 -scaling factor (mm/pixel) [1] := 3.1088 -matrix axis label [2] := y -!matrix size [2] := 97 -scaling factor (mm/pixel) [2] := 3.1088 -!slice thickness (pixels) := 3.375 -!END OF INTERFILE := diff --git a/recon_test_pack/lm_to_projdata.par b/recon_test_pack/lm_to_projdata.par index f2361ef1c..f5b1c291c 100644 --- a/recon_test_pack/lm_to_projdata.par +++ b/recon_test_pack/lm_to_projdata.par @@ -16,5 +16,6 @@ lm_to_projdata Parameters:= List event coordinates := 0 ; if you're short of RAM (i.e. a single projdata does not fit into memory), ; you can use this to process the list mode data in multiple passes. - num_segments_in_memory := -1 -End := +; num_segments_in_memory := 10 +End := + diff --git a/recon_test_pack/root_header.hroot b/recon_test_pack/root_header.hroot index d1056c2e0..390167ace 100644 --- a/recon_test_pack/root_header.hroot +++ b/recon_test_pack/root_header.hroot @@ -8,15 +8,13 @@ Inner ring diameter (cm) := 65.6 Average depth of interaction (cm) := 0.7 Distance between rings (cm) := 0.40625 Default bin size (cm) := 0.208626 -Maximum number of non-arc-corrected bins := 344 -Default number of arc-corrected bins := 344 +Maximum number of non-arc-corrected bins := 501 +Default number of arc-corrected bins := 501 View offset (degrees) := 0 -Number of TOF time bins := 411 -Size of timing bin (ps) := 10.00 +Number of TOF time bins := 5 +Size of timing bin (ps) := 820.00 Timing resolution (ps) := 400.0 -%TOF mashing factor:= 82 - GATE scanner type := GATE_Cylindrical_PET GATE_Cylindrical_PET Parameters := diff --git a/recon_test_pack/run_ML_norm_tests.sh b/recon_test_pack/run_ML_norm_tests.sh index bb4199f3b..73b3f342f 100755 --- a/recon_test_pack/run_ML_norm_tests.sh +++ b/recon_test_pack/run_ML_norm_tests.sh @@ -10,22 +10,14 @@ # Copyright (C) 2013, 2020, 2021 University College London # This file is part of STIR. # -# This file is free software; you can redistribute it and/or modify -# it under the terms of the GNU Lesser General Public License as published by -# the Free Software Foundation; either version 2.1 of the License, or -# (at your option) any later version. - -# This file is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU Lesser General Public License for more details. +# SPDX-License-Identifier: Apache-2.0 # # See STIR/LICENSE.txt for details # # Author Kris Thielemans # -echo This script should work with STIR version 5.0. If you have +echo This script should work with STIR version 5.2. If you have echo a later version, you might have to update your test pack. echo Please check the web site. echo diff --git a/recon_test_pack/run_root_GATE.sh b/recon_test_pack/run_root_GATE.sh index 0eb0f29c6..9b7edd5c5 100755 --- a/recon_test_pack/run_root_GATE.sh +++ b/recon_test_pack/run_root_GATE.sh @@ -11,7 +11,7 @@ # # Author Nikos Efthimiou, Kris Thielemans -echo This script should work with STIR version ">"3.0. If you have +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. echo Please check the web site. echo @@ -73,10 +73,7 @@ echo "GATE support detected!" fi -echo Executing tests on ROOT files generated by GATE simulations, with cylindrical PET scanners - - - +echo Executing tests on ROOT files generated by GATE simulations # first delete any files remaining from a previous run @@ -87,22 +84,12 @@ export INPUT_ROOT_FILE=test_PET_GATE.root export INPUT=root_header.hroot export TEMPLATE=template_for_ROOT_scanner.hs -# -# Get the number of events unlisted. -# -echo -echo ------------- Converting ROOT files to ProjData file ------------- -echo -echo Making ProjData for all events +echo "------------------------------" echo Running lm_to_projdata for all events export OUT_PROJDATA_FILE=my_proj_from_lm_all_events export EXCLUDE_SCATTERED=0 export EXCLUDE_RANDOM=0 -lm_to_projdata ./lm_to_projdata.par 2>"./my_root_log_all.log" -all_events=`awk -F ":" '/Number of prompts/ {print $2}' my_root_log_all.log` -echo "Number of prompts stored in this time period:" ${all_events} - log=lm_to_projdata_from_ROOT_all_events.log lm_to_projdata ./lm_to_projdata.par > ${log} 2>&1 if [ $? -ne 0 ]; then @@ -114,17 +101,11 @@ else all_events=$(cat ${log}|grep "Number of prompts stored in this time period" | grep -o -E '[0-9]+') echo Number of prompts stored in this time period: ${all_events} - echo fi -echo Reading all values from ROOT file -# -# Get the number of events directly from the ROOT file -# -echo -echo ------------- Reading values directly from ROOT file ------------- -echo -cat << EOF > my_root.input +echo Reading values directly from ROOT file +log=root_output_all_events.log +${ROOT} -b -l ${INPUT_ROOT_FILE} << EOF >& ${log} Coincidences->Draw(">>eventlist","","goff"); Int_t N = eventlist->GetN(); cout<& ${log} +${ROOT} -b -l ${INPUT_ROOT_FILE} << EOF >& ${log} Coincidences->Draw(">>eventlist","eventID1 == eventID2 && comptonPhantom1 == 0 && comptonPhantom2 == 0","goff"); Int_t N = eventlist->GetN(); cout<") + # target_link_libraries(IO PUBLIC PRIVATE "$") # So, we currently use an ugly work-around from # https://gitlab.kitware.com/cmake/cmake/-/issues/15415#note_334852 @@ -131,5 +131,5 @@ if (HAVE_JSON) endif() # currently needed for ParametricDensity (TODO get rid of this somehow?) -target_link_libraries(IO modelling_buildblock ) -target_link_libraries(IO listmode_buildblock) +target_link_libraries(IO PUBLIC modelling_buildblock ) +target_link_libraries(IO PUBLIC listmode_buildblock) diff --git a/src/IO/InterfileHeader.cxx b/src/IO/InterfileHeader.cxx index 7d3701762..484469478 100644 --- a/src/IO/InterfileHeader.cxx +++ b/src/IO/InterfileHeader.cxx @@ -306,12 +306,6 @@ bool InterfileHeader::post_processing() warning("Interfile error: no matrix size keywords present\n"); return true; } - if (matrix_size[matrix_size.size()-1].size()!=1) - { - warning("Interfile error: last dimension (%d) of 'matrix size' cannot be a list of numbers\n", - matrix_size[matrix_size.size()-1].size()); - return true; - } for (unsigned int dim=0; dim != matrix_size.size(); ++dim) { if (matrix_size[dim].size()==0) diff --git a/src/IO/InterfileHeaderSiemens.cxx b/src/IO/InterfileHeaderSiemens.cxx index 2e8f09b89..696be0d08 100644 --- a/src/IO/InterfileHeaderSiemens.cxx +++ b/src/IO/InterfileHeaderSiemens.cxx @@ -1,19 +1,16 @@ /* - Copyright (C) 2000 PARAPET partners - Copyright (C) 2000 - 2009-04-30, Hammersmith Imanet Ltd - Copyright (C) 2011-07-01 - 2012, Kris Thielemans - Copyright (C) 2013, 2016, 2018, 2020, 2023 University College London + Copyright (C) 2018, 2020, 2021, 2023 University College London Copyright (C) 2018 STFC This file is part of STIR. - SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license + SPDX-License-Identifier: Apache-2.0 See STIR/LICENSE.txt for details */ /*! \file \ingroup InterfileIO - \brief implementations for the stir::InterfileHeaderSiemens class + \brief implementations for the stir::InterfileHeaderSiemens classes \author Kris Thielemans \author PARAPET project @@ -29,6 +26,7 @@ #include "stir/IO/stir_ecat_common.h" #include #include +#include "stir/stream.h" #include "stir/warning.h" #include "stir/error.h" @@ -48,15 +46,17 @@ InterfileHeaderSiemens::InterfileHeaderSiemens() { // always PET exam_info_sptr->imaging_modality = ImagingModality::PT; - + + byte_order_values.clear(); byte_order_values.push_back("LITTLEENDIAN"); byte_order_values.push_back("BIGENDIAN"); - + + PET_data_type_values.clear(); PET_data_type_values.push_back("Emission"); PET_data_type_values.push_back("Transmission"); PET_data_type_values.push_back("Blank"); PET_data_type_values.push_back("AttenuationCorrection"); - PET_data_type_values.push_back("Normalisation"); + PET_data_type_values.push_back("Normalization"); PET_data_type_values.push_back("Image"); for (int patient_position_idx = 0; patient_position_idx <= PatientPosition::unknown_position; ++patient_position_idx) @@ -163,13 +163,43 @@ void InterfileHeaderSiemens::set_type_of_data() } +void +InterfileHeaderSiemens::ignore_Siemens_date_and_time_keys(const std::string& keyword) +{ + ignore_key(keyword + " date (yyyy:mm:dd)"); + ignore_key(keyword + " time (hh:mm:ss GMT+00:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT+01:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT+02:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT+03:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT+04:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT+05:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT+06:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT+07:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT+08:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT+09:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT+10:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT+11:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT-01:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT-02:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT-03:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT-04:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT-05:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT-06:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT-07:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT-08:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT-09:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT-10:00)"); + ignore_key(keyword + " time (hh:mm:ss GMT-11:00)"); +} + + /**********************************************************************/ InterfileRawDataHeaderSiemens::InterfileRawDataHeaderSiemens() : InterfileHeaderSiemens() { // first set to some crazy values - num_segments = -1; + num_segments = 0; num_rings = -1; maximum_ring_difference = -1; axial_compression = -1; @@ -208,13 +238,11 @@ InterfileRawDataHeaderSiemens::InterfileRawDataHeaderSiemens() ignore_key("%listmode header file"); ignore_key("%listmode data file"); ignore_key("%compressor version"); - ignore_key("%study date (yyyy"); - ignore_key("%study time (hh"); + ignore_Siemens_date_and_time_keys("%study"); ignore_key("isotope gamma halflife (sec)"); ignore_key("isotope branching factor"); ignore_key("radiopharmaceutical"); - ignore_key("%tracer injection date (yyyy"); - ignore_key("%tracer injection time (hh"); + ignore_Siemens_date_and_time_keys("%tracer injection"); ignore_key("relative time of tracer injection (sec)"); ignore_key("tracer activity at time of injection (bq)"); ignore_key("injected volume (ml)"); @@ -226,8 +254,7 @@ InterfileRawDataHeaderSiemens::InterfileRawDataHeaderSiemens() ignore_key("method of scatter correction"); ignore_key("%method of random correction"); ignore_key("%decay correction"); - ignore_key("%decay correction reference date (yyyy"); - ignore_key("%decay correction reference time (hh"); + ignore_Siemens_date_and_time_keys("%decay correction reference"); ignore_key("decay correction factor"); ignore_key("scatter fraction (%)"); ignore_key("scan data type description"); @@ -250,8 +277,8 @@ bool InterfileRawDataHeaderSiemens::post_processing() const std::string PET_data_type = standardise_interfile_keyword(PET_data_type_values[PET_data_type_index]); - if (PET_data_type != "emission" && PET_data_type != "transmission") - { error("Interfile error: expecting emission or transmission for 'PET data type'"); } + if (PET_data_type != "emission" && PET_data_type != "transmission" && PET_data_type != "normalization") + { error("Interfile error: expecting 'emission' or 'transmission' or 'normalization' for 'PET data type'"); } // handle scanner @@ -293,7 +320,7 @@ bool InterfileRawDataHeaderSiemens::post_processing() error("Interfile warning: 'number of segments' and length of 'segment table' are not consistent"); } segment_sequence = ecat::find_segment_sequence(*data_info_ptr); - //XXX check if order here and segment_table are consistent + //TODO check if order here and segment_table are consistent } // Set the bed position @@ -566,5 +593,115 @@ bool InterfileListmodeHeaderSiemens::post_processing() return false; } +InterfileNormHeaderSiemens::InterfileNormHeaderSiemens() + : InterfileRawDataHeaderSiemens() +{ + // some defaults + calib_factor = 1.F; + cross_calib_factor = 1.F; + num_buckets = 0; // should be set normally + num_components = 0; // should be set to 8 normally + axial_compression = 11; // should be set normally but seems to be this always + is_arccorrected = false; // norm data is never arc-corrected + + ignore_key("data description"); + ignore_Siemens_date_and_time_keys("%expiration"); + ignore_key("%raw normalization scans description"); + + // remove some standard keys, which Siemens has replaced with similar names + remove_key("matrix size"); + remove_key("matrix axis label"); + remove_key("scaling factor (mm/pixel)"); + + // keywords for the components + add_key("%number of normalization components", + KeyArgument::INT, (KeywordProcessor)&InterfileNormHeaderSiemens::read_num_components, &num_components); + add_vectorised_key("%matrix size", &matrix_size); + add_vectorised_key("%matrix axis label", &matrix_labels); + ignore_key("%matrix axis unit"); + ignore_key("%normalization component"); + ignore_key("%normalization components description"); + add_vectorised_key("data offset in bytes", &data_offset_each_dataset); + remove_key("number of dimensions"); + add_vectorised_key("number of dimensions", &number_of_dimensions); + ignore_key("%scale factor"); + + // other things + add_key("%number of buckets", &num_buckets); + + ignore_key("%global scanner calibration factor"); + add_key("%scanner quantification factor (Bq*s/ECAT counts)",& calib_factor); + add_key("%cross calibration factor",& cross_calib_factor); + ignore_Siemens_date_and_time_keys("%calibration"); + + // isotope things are vectorised in norm files and not in other raw data, so we could + // fix that, but as we are not interested in it anyway (tends to be Ge-68), let's just ignore it. + remove_key("isotope name"); + ignore_key("isotope name"); + + ignore_key("%number of normalization scans"); + ignore_key("%normalization scan"); + remove_key("image duration (sec)"); + ignore_key("image duration (sec)"); +#if 0 + // change to vectorised key + // would need to set image_durations length from "number of normalization scans" + add_vectorised_key("image duration (sec)", &image_durations); +#endif + ignore_key("%data format"); + ignore_key("%data set description"); + ignore_key("total number of data sets"); + ignore_key("%data set"); +} + +void InterfileNormHeaderSiemens::read_num_components() +{ + set_variable(); + + matrix_labels.resize(num_components); + matrix_size.resize(num_components); + // matrix_axis_units.resize(num_components); + // matrix_axis_units.resize(num_components); + // pixel_sizes.resize(num_components); + // normalization_components.resize(num_components); + data_offset_each_dataset.resize(num_components); + number_of_dimensions.resize(num_components); +} + +bool InterfileNormHeaderSiemens::post_processing() +{ + if (matrix_size.size() < 4) + error("Error parsing ECAT8 norm file header: '%number of normalization components='" + + std::to_string(matrix_size.size())+ + "' but should be at least 4"); + // %normalization component [1]:=geometric effects + if (matrix_size[0].size() != 2) + error("Error parsing ECAT8 norm file header: '%matrix size[1]' should have length 2"); + // %normalization component [3]:=crystal efficiencies + if (matrix_size[2].size() != 2) + error("Error parsing ECAT8 norm file header: '%matrix size[3]' should have length 2"); + + // TODO should do far more checks... + + // remove trailing \r (or other white space) occuring in mMR norm files (they sometimes have \r\r at end of line) + std::string s=exam_info_sptr->originating_system; + s.erase( std::remove_if( s.begin(), s.end(), isspace ), s.end() ); + exam_info_sptr->originating_system=s; + s=data_file_name; + s.erase( std::remove_if( s.begin(), s.end(), isspace ), s.end() ); + data_file_name=s; + + // norm headers don't seem to have "number of views". We need to get it from elsewhere... + // The crystal efficiencies have as first dimension the number of crystals, so let's use that. + const int num_detectors_per_ring = matrix_size[2][0]; + this->num_views = num_detectors_per_ring/2; + // find num_bins from geometric effects + this->num_bins = matrix_size[0][0]; + + if (InterfileRawDataHeaderSiemens::post_processing() == true) + return true; + + return false; +} END_NAMESPACE_STIR diff --git a/src/IO/interfile.cxx b/src/IO/interfile.cxx index 0b129cdab..dd77e45d1 100644 --- a/src/IO/interfile.cxx +++ b/src/IO/interfile.cxx @@ -56,9 +56,8 @@ #include #include "stir/ProjDataInfoBlocksOnCylindricalNoArcCorr.h" #include "stir/ProjDataInfoGenericNoArcCorr.h" +#include "stir/ProjDataInfoSubsetByView.h" - -#ifndef STIR_NO_NAMESPACES using std::cerr; using std::endl; using std::ofstream; @@ -70,7 +69,6 @@ using std::vector; using std::istream; using std::string; using std::ios; -#endif START_NAMESPACE_STIR @@ -488,7 +486,7 @@ static void write_interfile_image_data_descriptions(std::ostream& output_header, output_header << "number of image data types := " << data_type_descriptions.size() << '\n'; output_header << "index nesting level := {data type}\n"; - for (int i=0; i file_offsets(image.get_num_params()); VectorWithOffset scaling_factors(image.get_num_params()); - for (int i=1; i<=image.get_num_params(); i++) { + for (int i=1; i<=static_cast(image.get_num_params()); i++) { float scale_to_use = scale; file_offsets[i-1] = output_data.tellp(); write_data(output_data, image.construct_single_density(i), output_type, scale_to_use, @@ -958,7 +956,7 @@ write_basic_interfile(const string& filename, VectorWithOffset file_offsets(image.get_num_time_frames()); VectorWithOffset scaling_factors(image.get_num_time_frames()); - for (int i=1; i<=image.get_num_time_frames(); i++) + for (int i=1; i<=static_cast(image.get_num_time_frames()); i++) { float scale_to_use = scale; file_offsets[i-1] = output_data.tellp(); @@ -1535,8 +1533,14 @@ write_basic_interfile_PDFS_header(const string& header_file_name, << proj_data_info_sptr->get_sampling_in_s(Bin(0,0,0,0))/10. << endl; } // end generic scanner - - else error("write_basic_interfile_PDFS_header: Error casting the projdata to one of its geometries: Cylindrical/BlocksOnCylindrical/Genreic"); + else if (!dynamic_pointer_cast(pdfs.get_proj_data_info_sptr())) + { + error("write_basic_interfile_PDFS_header: cannot write subset data yet. Sorry"); + } + else + { + error("write_basic_interfile_PDFS_header: Error casting the projdata to one of its geometries: Cylindrical/BlocksOnCylindrical/Generic"); + } } } diff --git a/src/Shape_buildblock/CMakeLists.txt b/src/Shape_buildblock/CMakeLists.txt index 61e3a5123..05b813601 100644 --- a/src/Shape_buildblock/CMakeLists.txt +++ b/src/Shape_buildblock/CMakeLists.txt @@ -18,4 +18,4 @@ set(${dir_LIB_SOURCES} include(stir_lib_target) -target_link_libraries(Shape_buildblock buildblock IO numerics_buildblock ) +target_link_libraries(Shape_buildblock PUBLIC buildblock IO numerics_buildblock ) diff --git a/src/analytic/FBP2D/CMakeLists.txt b/src/analytic/FBP2D/CMakeLists.txt index d72dce63e..533520f44 100644 --- a/src/analytic/FBP2D/CMakeLists.txt +++ b/src/analytic/FBP2D/CMakeLists.txt @@ -13,7 +13,7 @@ set(${dir_LIB_SOURCES} #$(dir)_REGISTRY_SOURCES:= include(stir_lib_target) -target_link_libraries(analytic_FBP2D recon_buildblock IO ) +target_link_libraries(analytic_FBP2D PUBLIC recon_buildblock IO ) set (dir_EXE_SOURCES ${dir}_EXE_SOURCES) diff --git a/src/analytic/FBP3DRP/CMakeLists.txt b/src/analytic/FBP3DRP/CMakeLists.txt index 892b14c44..860bb5156 100644 --- a/src/analytic/FBP3DRP/CMakeLists.txt +++ b/src/analytic/FBP3DRP/CMakeLists.txt @@ -13,7 +13,7 @@ set(${dir_LIB_SOURCES} #$(dir)_REGISTRY_SOURCES:= include(stir_lib_target) -target_link_libraries(analytic_FBP3DRP analytic_FBP2D recon_buildblock ) +target_link_libraries(analytic_FBP3DRP PUBLIC analytic_FBP2D recon_buildblock ) set (dir_EXE_SOURCES ${dir}_EXE_SOURCES) diff --git a/src/buildblock/CMakeLists.txt b/src/buildblock/CMakeLists.txt index 7842ba777..1d90c88c6 100644 --- a/src/buildblock/CMakeLists.txt +++ b/src/buildblock/CMakeLists.txt @@ -130,17 +130,17 @@ endif() # TODO Remove but currently needed for ProjData.cxx, DynamicDisc*cxx, TimeFrameDef if (LLN_FOUND) - target_link_libraries(buildblock ${LLN_LIBRARIES}) + target_link_libraries(buildblock PUBLIC ${LLN_LIBRARIES}) endif() if (RDF_FOUND) # TODO cannot do this as it creates circular dependencies - # target_link_libraries(buildblock local_IO_GE) + # target_link_libraries(buildblock PUBLIC local_IO_GE) endif() # TODO currently needed as filters need fourier -#target_link_libraries(buildblock numerics_buildblock) +#target_link_libraries(buildblock PUBLIC numerics_buildblock) if (STIR_OPENMP) - target_link_libraries(buildblock ${OpenMP_EXE_LINKER_FLAGS}) + target_link_libraries(buildblock PUBLIC ${OpenMP_EXE_LINKER_FLAGS}) endif() diff --git a/src/buildblock/KeyParser.cxx b/src/buildblock/KeyParser.cxx index da19402a7..ab609755a 100644 --- a/src/buildblock/KeyParser.cxx +++ b/src/buildblock/KeyParser.cxx @@ -309,8 +309,12 @@ string KeyParser::get_keyword(const string& line) const { // keyword stops at either := or an index [] - // TODO should check that = follows : to allow keywords with colons in there - const string::size_type eok = line.find_first_of(":[",0); + auto eok = line.find_first_of(":[",0); + // check that = follows : to allow keywords containing colons + while (eok != string::npos && line[eok] == ':' && eok+1 < line.size() && line[eok+1] != '=') + { + eok = line.find_first_of(":[", eok+1); + } return line.substr(0,eok); } diff --git a/src/buildblock/ProjData.cxx b/src/buildblock/ProjData.cxx index c021806ca..c8582bab0 100644 --- a/src/buildblock/ProjData.cxx +++ b/src/buildblock/ProjData.cxx @@ -3,7 +3,7 @@ Copyright (C) 2000 - 2010-10-15, Hammersmith Imanet Ltd Copyright (C) 2011-07-01 -2013, Kris Thielemans Copyright (C) 2016, University of Hull - Copyright (C) 2015, 2020, 2022 University College London + Copyright (C) 2015, 2020, 2022, 2023 University College London Copyright (C) 2021-2022, Commonwealth Scientific and Industrial Research Organisation Copyright (C) 2021, Rutherford Appleton Laboratory STFC This file is part of STIR. @@ -60,21 +60,18 @@ #include "stir/IO/GEHDF5Wrapper.h" #endif #include "stir/IO/stir_ecat7.h" -#include "stir/ViewSegmentNumbers.h" +#include "stir/ViewgramIndices.h" #include "stir/is_null_ptr.h" #include #include #include -#include "stir/warning.h" #include "stir/error.h" -#ifndef STIR_NO_NAMESPACES using std::istream; using std::fstream; using std::ios; using std::string; using std::vector; -#endif START_NAMESPACE_STIR @@ -130,8 +127,7 @@ read_from_file(const string& filename, #ifndef STIR_USE_GE_IO { #ifndef NDEBUG - warning("ProjData::read_from_file trying to read %s as GE Advance file", - filename.c_str()); + info("ProjData::read_from_file trying to read " + filename + " as GE Advance file", 3); #endif return shared_ptr( new ProjDataGEAdvance(input) ); } @@ -139,8 +135,7 @@ read_from_file(const string& filename, #else // use VOLPET { #ifndef NDEBUG - warning("ProjData::read_from_file trying to read %s as GE VOLPET file", - filename.c_str()); + info("ProjData::read_from_file trying to read " + filename + " as GE VOLPET file", 3); #endif delete input;// TODO no longer use pointer after getting rid of ProjDataGEAdvance return shared_ptr( new GE_IO::ProjDataVOLPET(filename, openmode) ); @@ -155,8 +150,7 @@ read_from_file(const string& filename, if (GE_IO::is_IE_signature(signature)) { #ifndef NDEBUG - warning("ProjData::read_from_file trying to read %s as GE IE file", - filename.c_str()); + info("ProjData::read_from_file trying to read " + filename + " as GE IE file", 3); #endif return shared_ptr( new GE_IO::ProjDataIE(filename) ); } @@ -168,22 +162,21 @@ read_from_file(const string& filename, if (strncmp(signature, "MATRIX", 6) == 0) { #ifndef NDEBUG - warning("ProjData::read_from_file trying to read %s as ECAT7", filename.c_str()); + info("ProjData::read_from_file trying to read " + filename + " as ECAT7", 3); #endif USING_NAMESPACE_ECAT; USING_NAMESPACE_ECAT7; if (is_ECAT7_emission_file(actual_filename) || is_ECAT7_attenuation_file(actual_filename)) { - warning("\nReading frame 1, gate 1, data 0, bed 0 from file %s", - actual_filename.c_str()); + info("Reading frame 1, gate 1, data 0, bed 0 from file " + actual_filename, 3); shared_ptr proj_data_sptr(ECAT7_to_PDFS(filename, /*frame_num, gate_num, data_num, bed_num*/1,1,0,0)); return proj_data_sptr; } else { if (is_ECAT7_file(actual_filename)) - warning("ProjData::read_from_file ECAT7 file %s is of unsupported file type", actual_filename.c_str()); + error("ProjData::read_from_file ECAT7 file " + actual_filename + " is of unsupported file type"); } } #endif // HAVE_LLN_MATRIX @@ -192,7 +185,7 @@ read_from_file(const string& filename, if (is_interfile_signature(signature)) { #ifndef NDEBUG - warning("ProjData::read_from_file trying to read %s as Interfile", filename.c_str()); + info("ProjData::read_from_file trying to read " + filename + " as Interfile", 3); #endif shared_ptr ptr(read_interfile_PDFS(filename, openmode)); if (!is_null_ptr(ptr)) @@ -204,7 +197,7 @@ read_from_file(const string& filename, if (GE_IO::is_RDF_file(actual_filename)) { #ifndef NDEBUG - warning("ProjData::read_from_file trying to read %s as RDF", filename.c_str()); + info("ProjData::read_from_file trying to read " + filename + " as RDF", 3); #endif shared_ptr ptr(new GE_IO::ProjDataRDF(filename)); if (!is_null_ptr(ptr)) @@ -216,7 +209,7 @@ read_from_file(const string& filename, if (GE::RDF_HDF5::GEHDF5Wrapper::check_GE_signature(actual_filename)) { #ifndef NDEBUG - warning("ProjData::read_from_file trying to read %s as GE HDF5", filename.c_str()); + info("ProjData::read_from_file trying to read " + filename + " as GE HDF5", 3); #endif shared_ptr ptr(new GE::RDF_HDF5::ProjDataGEHDF5(filename)); if (!is_null_ptr(ptr)) @@ -224,9 +217,8 @@ read_from_file(const string& filename, } #endif // GE HDF5 - error("\nProjData::read_from_file could not read projection data %s.\n" - "Unsupported file format? Aborting.", - filename.c_str()); + error("ProjData::read_from_file could not read projection data " + filename + ".\n" + "Unsupported file format? Aborting."); // need to return something to satisfy the compiler, but we never get here shared_ptr null_ptr; return null_ptr; @@ -263,7 +255,19 @@ ProjData::get_subset(const std::vector& views) const } -Viewgram +Viewgram +ProjData::get_empty_viewgram(const ViewgramIndices& ind) const +{ + return proj_data_info_sptr->get_empty_viewgram(ind); +} + +Sinogram +ProjData::get_empty_sinogram(const SinogramIndices& ind) const +{ + return proj_data_info_sptr->get_empty_sinogram(ind); +} + +Viewgram ProjData::get_empty_viewgram(const int view_num, const int segment_num, const bool make_num_tangential_poss_odd, const int timing_pos) const @@ -301,12 +305,24 @@ ProjData::get_empty_segment_by_view(const int segment_num, proj_data_info_sptr->get_empty_segment_by_view(segment_num, make_num_tangential_poss_odd, timing_pos); } -RelatedViewgrams -ProjData::get_empty_related_viewgrams(const ViewSegmentNumbers& view_segmnet_num, - //const int view_num, const int segment_num, - const shared_ptr& symmetries_used, - const bool make_num_tangential_poss_odd, - const int timing_pos) const +SegmentBySinogram +ProjData::get_empty_segment_by_sinogram(const SegmentIndices& ind) const +{ + return proj_data_info_sptr->get_empty_segment_by_sinogram(ind); +} + +SegmentByView +ProjData::get_empty_segment_by_view(const SegmentIndices& ind) const +{ + return proj_data_info_sptr->get_empty_segment_by_view(ind); +} + +RelatedViewgrams +ProjData::get_empty_related_viewgrams(const ViewgramIndices& view_segmnet_num, + const shared_ptr& symmetries_used, + const bool make_num_tangential_poss_odd, + const int timing_pos) const + { return proj_data_info_sptr->get_empty_related_viewgrams(view_segmnet_num, symmetries_used, make_num_tangential_poss_odd, timing_pos); @@ -314,16 +330,15 @@ ProjData::get_empty_related_viewgrams(const ViewSegmentNumbers& view_segmnet_num RelatedViewgrams -ProjData::get_related_viewgrams(const ViewSegmentNumbers& view_segment_num, - //const int view_num, const int segment_num, - const shared_ptr& symmetries_used, - const bool make_num_bins_odd, - const int timing_pos) const +ProjData::get_related_viewgrams(const ViewgramIndices& viewgram_indices, + const shared_ptr& symmetries_used, + const bool make_num_tangential_poss_odd, + const int timing_pos) const { vector pairs; symmetries_used->get_related_view_segment_numbers( pairs, - ViewSegmentNumbers(view_segment_num.view_num(),view_segment_num.segment_num()) + viewgram_indices ); vector > viewgrams; @@ -332,8 +347,9 @@ ProjData::get_related_viewgrams(const ViewSegmentNumbers& view_segment_num, for (unsigned int i=0; iget_viewgram(pairs[i])); } return RelatedViewgrams(viewgrams, symmetries_used); diff --git a/src/buildblock/ProjDataInfo.cxx b/src/buildblock/ProjDataInfo.cxx index f032536f0..9f6ab2542 100644 --- a/src/buildblock/ProjDataInfo.cxx +++ b/src/buildblock/ProjDataInfo.cxx @@ -5,7 +5,7 @@ Copyright (C) 2000 - 2009-05-13, Hammersmith Imanet Ltd Copyright (C) 2011-07-01 - 2011, Kris Thielemans Copyright (C) 2018, University of Leeds - Copyright (C) 2018, 2020-2022 University College London + Copyright (C) 2018, 2020-2023 University College London Copyright (C) 2016-2019, University of Hull This file is part of STIR. @@ -394,6 +394,15 @@ ProjDataInfo::get_empty_viewgram(const int view_num, return v; } +Viewgram +ProjDataInfo::get_empty_viewgram(const ViewgramIndices& ind) const +{ + // we can't access the shared ptr, so we have to clone 'this'. + shared_ptr proj_data_info_sptr(this->clone()); + Viewgram v(proj_data_info_sptr, ind); + return v; +} + Sinogram ProjDataInfo::get_empty_sinogram(const int axial_pos_num, const int segment_num, @@ -411,6 +420,15 @@ ProjDataInfo::get_empty_sinogram(const int axial_pos_num, const int segment_num, return s; } +Sinogram +ProjDataInfo::get_empty_sinogram(const SinogramIndices& ind) const +{ + // we can't access the shared ptr, so we have to clone 'this'. + shared_ptr proj_data_info_sptr(this->clone()); + Sinogram s(proj_data_info_sptr, ind); + return s; +} + SegmentBySinogram ProjDataInfo::get_empty_segment_by_sinogram(const int segment_num, const bool make_num_tangential_poss_odd, @@ -430,6 +448,14 @@ ProjDataInfo::get_empty_segment_by_sinogram(const int segment_num, return s; } +SegmentBySinogram +ProjDataInfo::get_empty_segment_by_sinogram(const SegmentIndices& ind) const +{ + // we can't access the shared ptr, so we have to clone 'this'. + shared_ptr proj_data_info_sptr(this->clone()); + SegmentBySinogram s(proj_data_info_sptr, ind); + return s; +} SegmentByView ProjDataInfo::get_empty_segment_by_view(const int segment_num, @@ -450,29 +476,35 @@ ProjDataInfo::get_empty_segment_by_view(const int segment_num, return s; } +SegmentByView +ProjDataInfo::get_empty_segment_by_view(const SegmentIndices& ind) const +{ + // we can't access the shared ptr, so we have to clone 'this'. + shared_ptr proj_data_info_sptr(this->clone()); + SegmentByView s(proj_data_info_sptr, ind); + return s; +} + RelatedViewgrams -ProjDataInfo::get_empty_related_viewgrams(const ViewSegmentNumbers& view_segment_num, - //const int view_num, const int segment_num, +ProjDataInfo::get_empty_related_viewgrams(const ViewgramIndices& viewgram_indices, const shared_ptr& symmetries_used, const bool make_num_tangential_poss_odd, const int timing_pos_num) const { + if (make_num_tangential_poss_odd) + error("make_num_tangential_poss_odd is no longer supported"); vector pairs; - symmetries_used->get_related_view_segment_numbers( - pairs, - ViewSegmentNumbers(view_segment_num.view_num(),view_segment_num.segment_num()) - ); + symmetries_used->get_related_view_segment_numbers(pairs, viewgram_indices); vector > viewgrams; viewgrams.reserve(pairs.size()); for (unsigned int i=0; i(viewgrams, symmetries_used); @@ -750,8 +782,9 @@ blindly_equals(const root_type * const that) const (get_max_view_num()==proj.get_max_view_num()) && (get_min_tangential_pos_num() ==proj.get_min_tangential_pos_num())&& (get_max_tangential_pos_num() ==proj.get_max_tangential_pos_num())&& - (proj.is_tof_data() && is_tof_data() ? (get_min_tof_pos_num() == proj.get_min_tof_pos_num()) && - (get_max_tof_pos_num() == proj.get_max_tof_pos_num()) : true) && + (get_tof_mash_factor() == proj.get_tof_mash_factor()) && + (get_min_tof_pos_num() == proj.get_min_tof_pos_num()) && + (get_max_tof_pos_num() == proj.get_max_tof_pos_num()) && equal(min_axial_pos_per_seg.begin(), min_axial_pos_per_seg.end(), proj.min_axial_pos_per_seg.begin())&& equal(max_axial_pos_per_seg.begin(), max_axial_pos_per_seg.end(), proj.max_axial_pos_per_seg.begin())&& (*get_scanner_ptr()== *(proj.get_scanner_ptr()))&& @@ -782,7 +815,7 @@ operator !=(const root_type& that) const \c true only if the types are the same, they are equal, or the range for the TOF, segments, axial and tangential positions is at least as large. - \warning Currently view ranges have to be identical. + \warning Currently view and TOF ranges have to be identical. */ bool ProjDataInfo:: @@ -796,13 +829,15 @@ operator>=(const ProjDataInfo& proj_data_info) const if (larger_proj_data_info == proj_data_info) return true; + if (proj_data_info.get_tof_mash_factor() != larger_proj_data_info.get_tof_mash_factor()) + return false; + if (proj_data_info.get_max_segment_num() > larger_proj_data_info.get_max_segment_num() || proj_data_info.get_min_segment_num() < larger_proj_data_info.get_min_segment_num() || proj_data_info.get_max_tangential_pos_num() > larger_proj_data_info.get_max_tangential_pos_num() || proj_data_info.get_min_tangential_pos_num() < larger_proj_data_info.get_min_tangential_pos_num() || - ((proj_data_info.get_min_tof_pos_num() < larger_proj_data_info.get_min_tof_pos_num() || - proj_data_info.get_max_tof_pos_num() > larger_proj_data_info.get_max_tof_pos_num()) && - (proj_data_info.is_tof_data() && larger_proj_data_info.is_tof_data()))) + proj_data_info.get_min_tof_pos_num() < larger_proj_data_info.get_min_tof_pos_num() || + proj_data_info.get_max_tof_pos_num() > larger_proj_data_info.get_max_tof_pos_num()) return false; for (int segment_num=proj_data_info.get_min_segment_num(); diff --git a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx index 1c101e24b..f3cbe6e77 100644 --- a/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx +++ b/src/buildblock/ProjDataInfoCylindricalNoArcCorr.cxx @@ -32,22 +32,15 @@ #include "stir/round.h" #include #include "stir/error.h" - -#ifdef BOOST_NO_STRINGSTREAM -#include -#else #include -#endif #include -#ifndef STIR_NO_NAMESPACES using std::endl; using std::ends; using std::string; using std::pair; using std::vector; -#endif START_NAMESPACE_STIR ProjDataInfoCylindricalNoArcCorr:: @@ -55,26 +48,28 @@ ProjDataInfoCylindricalNoArcCorr() {} ProjDataInfoCylindricalNoArcCorr:: -ProjDataInfoCylindricalNoArcCorr(const shared_ptr scanner_ptr, +ProjDataInfoCylindricalNoArcCorr(const shared_ptr scanner_sptr, const float ring_radius_v, const float angular_increment_v, const VectorWithOffset& num_axial_pos_per_segment, const VectorWithOffset& min_ring_diff_v, const VectorWithOffset& max_ring_diff_v, const int num_views,const int num_tangential_poss, const int tof_mash_factor) -: ProjDataInfoCylindrical(scanner_ptr, +: ProjDataInfoCylindrical(scanner_sptr, num_axial_pos_per_segment, min_ring_diff_v, max_ring_diff_v, num_views, num_tangential_poss), ring_radius(ring_radius_v), angular_increment(angular_increment_v) { - if (num_tangential_poss > scanner_ptr->get_max_num_non_arccorrected_bins()) - error("Configured tangential positions exceed the maximum number of non arc-corrected bins set for the scanner."); - + if (!scanner_sptr) + error("ProjDataInfoCylindricalNoArcCorr: first argument (scanner_ptr) is zero"); + if (num_tangential_poss > scanner_sptr->get_max_num_non_arccorrected_bins()) + error("ProjDataInfoCylindricalNoArcCorr: number of tangential positions exceeds the maximum number of non arc-corrected bins set for the scanner."); + uncompressed_view_tangpos_to_det1det2_initialised = false; det1det2_to_uncompressed_view_tangpos_initialised = false; - if(scanner_ptr->is_tof_ready()) + if(scanner_sptr->is_tof_ready()) set_tof_mash_factor(tof_mash_factor); #ifdef STIR_OPENMP_SAFE_BUT_SLOW this->initialise_uncompressed_view_tangpos_to_det1det2(); @@ -83,30 +78,20 @@ ProjDataInfoCylindricalNoArcCorr(const shared_ptr scanner_ptr, } ProjDataInfoCylindricalNoArcCorr:: -ProjDataInfoCylindricalNoArcCorr(const shared_ptr scanner_ptr, +ProjDataInfoCylindricalNoArcCorr(const shared_ptr scanner_sptr, const VectorWithOffset& num_axial_pos_per_segment, const VectorWithOffset& min_ring_diff_v, const VectorWithOffset& max_ring_diff_v, const int num_views, const int num_tangential_poss, const int tof_mash_factor) -: ProjDataInfoCylindrical(scanner_ptr, + : ProjDataInfoCylindricalNoArcCorr(scanner_sptr, + scanner_sptr ? scanner_sptr->get_effective_ring_radius() : 0.F, // avoid segfault if scanner_sptr==0 + scanner_sptr ? static_cast(_PI/scanner_sptr->get_num_detectors_per_ring()) : 0.F, num_axial_pos_per_segment, min_ring_diff_v, max_ring_diff_v, - num_views, num_tangential_poss) -{ - assert(!is_null_ptr(scanner_ptr)); - ring_radius = scanner_ptr->get_effective_ring_radius(); - angular_increment = static_cast(_PI/scanner_ptr->get_num_detectors_per_ring()); - uncompressed_view_tangpos_to_det1det2_initialised = false; - det1det2_to_uncompressed_view_tangpos_initialised = false; - - if(scanner_ptr->is_tof_ready()) - set_tof_mash_factor(tof_mash_factor); -#ifdef STIR_OPENMP_SAFE_BUT_SLOW - this->initialise_uncompressed_view_tangpos_to_det1det2(); - this->initialise_det1det2_to_uncompressed_view_tangpos(); -#endif -} + num_views, num_tangential_poss, + tof_mash_factor) +{} diff --git a/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx b/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx index 6ea435619..e900897d2 100644 --- a/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx +++ b/src/buildblock/ProjDataInfoGenericNoArcCorr.cxx @@ -52,20 +52,21 @@ ProjDataInfoGenericNoArcCorr() {} ProjDataInfoGenericNoArcCorr:: -ProjDataInfoGenericNoArcCorr(const shared_ptr scanner_ptr, +ProjDataInfoGenericNoArcCorr(const shared_ptr scanner_sptr, const VectorWithOffset& num_axial_pos_per_segment, const VectorWithOffset& min_ring_diff_v, const VectorWithOffset& max_ring_diff_v, const int num_views,const int num_tangential_poss) -: ProjDataInfoGeneric(scanner_ptr, +: ProjDataInfoGeneric(scanner_sptr, num_axial_pos_per_segment, min_ring_diff_v, max_ring_diff_v, num_views, num_tangential_poss) { - if (num_tangential_poss > scanner_ptr->get_max_num_non_arccorrected_bins()) - error("Configured tangential positions exceed the maximum number of non arc-corrected bins set for the scanner."); + if (!scanner_sptr) + error("ProjDataInfoGenericNoArcCorr: first argument (scanner_ptr) is zero"); + if (num_tangential_poss > scanner_sptr->get_max_num_non_arccorrected_bins()) + error("ProjDataInfoGenericNoArcCorr: number of tangential positions exceeds the maximum number of non arc-corrected bins set for the scanner."); - assert(!is_null_ptr(scanner_ptr)); uncompressed_view_tangpos_to_det1det2_initialised = false; det1det2_to_uncompressed_view_tangpos_initialised = false; #ifdef STIR_OPENMP_SAFE_BUT_SLOW diff --git a/src/buildblock/SegmentBySinogram.cxx b/src/buildblock/SegmentBySinogram.cxx index e9d0f58db..952d6a4a1 100644 --- a/src/buildblock/SegmentBySinogram.cxx +++ b/src/buildblock/SegmentBySinogram.cxx @@ -34,42 +34,57 @@ START_NAMESPACE_STIR template SegmentBySinogram :: SegmentBySinogram(const Array<3,elemT>& v, - const shared_ptr& pdi_ptr, - const int segment_num, - const int timing_pos_num) + const shared_ptr& pdi_ptr, + const SegmentIndices& ind) : - Segment(pdi_ptr, segment_num, timing_pos_num), + Segment(pdi_ptr, ind), Array<3,elemT>(v) { assert( get_min_view_num() == pdi_ptr->get_min_view_num()); assert( get_max_view_num() == pdi_ptr->get_max_view_num()); - assert( get_min_axial_pos_num() == pdi_ptr->get_min_axial_pos_num(segment_num)); - assert( get_max_axial_pos_num() == pdi_ptr->get_max_axial_pos_num(segment_num)); + assert( get_min_axial_pos_num() == pdi_ptr->get_min_axial_pos_num(ind.segment_num())); + assert( get_max_axial_pos_num() == pdi_ptr->get_max_axial_pos_num(ind.segment_num())); assert( get_min_tangential_pos_num() == pdi_ptr->get_min_tangential_pos_num()); assert( get_max_tangential_pos_num() == pdi_ptr->get_max_tangential_pos_num()); } template SegmentBySinogram :: -SegmentBySinogram(const shared_ptr &pdi_ptr, - const int segment_num, - const int timing_pos_num) +SegmentBySinogram(const shared_ptr& pdi_ptr, + const SegmentIndices& ind) : - Segment(pdi_ptr, segment_num, timing_pos_num), - Array<3,elemT>(IndexRange3D(pdi_ptr->get_min_axial_pos_num(segment_num), - pdi_ptr->get_max_axial_pos_num(segment_num), + Segment(pdi_ptr, ind), + Array<3,elemT>(IndexRange3D(pdi_ptr->get_min_axial_pos_num(ind.segment_num()), + pdi_ptr->get_max_axial_pos_num(ind.segment_num()), pdi_ptr->get_min_view_num(), pdi_ptr->get_max_view_num(), pdi_ptr->get_min_tangential_pos_num(), pdi_ptr->get_max_tangential_pos_num())) {} +template +SegmentBySinogram:: +SegmentBySinogram(const Array<3,elemT>& v, + const shared_ptr& pdi_sptr, + int segment_num, int timing_pos_num) + : + SegmentBySinogram(v, pdi_sptr, SegmentIndices(segment_num, timing_pos_num)) +{} + +template +SegmentBySinogram:: +SegmentBySinogram(const shared_ptr& pdi_sptr, + const int segment_num, const int t_num) + : + SegmentBySinogram(pdi_sptr, SegmentIndices(segment_num, t_num)) +{} + template SegmentBySinogram:: SegmentBySinogram(const SegmentByView& s_v ) : Segment(s_v.get_proj_data_info_sptr()->create_shared_clone(), - s_v.get_segment_num(), s_v.get_timing_pos_num()), + s_v.get_segment_indices()), Array<3,elemT> (IndexRange3D (s_v.get_min_axial_pos_num(), s_v.get_max_axial_pos_num(), s_v.get_min_view_num(), s_v.get_max_view_num(), s_v.get_min_tangential_pos_num(), s_v.get_max_tangential_pos_num())) diff --git a/src/buildblock/SegmentByView.cxx b/src/buildblock/SegmentByView.cxx index 402f64cc8..7e5dd443e 100644 --- a/src/buildblock/SegmentByView.cxx +++ b/src/buildblock/SegmentByView.cxx @@ -32,17 +32,16 @@ START_NAMESPACE_STIR template SegmentByView:: SegmentByView(const Array<3,elemT>& v, - const shared_ptr& pdi_ptr, - const int segment_num, - const int timing_pos_num) + const shared_ptr& pdi_ptr, + const SegmentIndices& ind) : - Segment(pdi_ptr, segment_num, timing_pos_num), + Segment(pdi_ptr, ind), Array<3,elemT>(v) { assert( get_min_view_num() == pdi_ptr->get_min_view_num()); assert( get_max_view_num() == pdi_ptr->get_max_view_num()); - assert( get_min_axial_pos_num() == pdi_ptr->get_min_axial_pos_num(segment_num)); - assert( get_max_axial_pos_num() == pdi_ptr->get_max_axial_pos_num(segment_num)); + assert( get_min_axial_pos_num() == pdi_ptr->get_min_axial_pos_num(ind.segment_num())); + assert( get_max_axial_pos_num() == pdi_ptr->get_max_axial_pos_num(ind.segment_num())); assert( get_min_tangential_pos_num() == pdi_ptr->get_min_tangential_pos_num()); assert( get_max_tangential_pos_num() == pdi_ptr->get_max_tangential_pos_num()); @@ -51,23 +50,37 @@ SegmentByView(const Array<3,elemT>& v, template SegmentByView:: SegmentByView(const shared_ptr& pdi_ptr, - const int segment_num, - const int timing_pos_num) + const SegmentIndices& ind) : - Segment(pdi_ptr, segment_num, timing_pos_num), + Segment(pdi_ptr, ind), Array<3,elemT>(IndexRange3D(pdi_ptr->get_min_view_num(), pdi_ptr->get_max_view_num(), - pdi_ptr->get_min_axial_pos_num(segment_num), - pdi_ptr->get_max_axial_pos_num(segment_num), + pdi_ptr->get_min_axial_pos_num(ind.segment_num()), + pdi_ptr->get_max_axial_pos_num(ind.segment_num()), pdi_ptr->get_min_tangential_pos_num(), pdi_ptr->get_max_tangential_pos_num())) {} +template +SegmentByView:: +SegmentByView(const Array<3,elemT>& v, + const shared_ptr& pdi_sptr, + const int segment_num, const int t_num) + : + SegmentByView(v, pdi_sptr, SegmentIndices(segment_num, t_num)) +{} + +template +SegmentByView:: +SegmentByView(const shared_ptr& pdi_sptr, + const int segment_num, const int t_num) + : SegmentByView(pdi_sptr, SegmentIndices(segment_num, t_num)) +{} + template SegmentByView::SegmentByView(const SegmentBySinogram& s_s) : Segment(s_s.get_proj_data_info_sptr()->create_shared_clone(), - s_s.get_segment_num(),s_s.get_timing_pos_num()), - + s_s.get_segment_indices()), Array<3,elemT> (IndexRange3D(s_s.get_min_view_num(),s_s.get_max_view_num(), s_s.get_min_axial_pos_num(),s_s.get_max_axial_pos_num(), s_s.get_min_tangential_pos_num(), s_s.get_max_tangential_pos_num())) diff --git a/src/buildblock/extend_projdata.cxx b/src/buildblock/extend_projdata.cxx index 1258abd8e..774a62e96 100644 --- a/src/buildblock/extend_projdata.cxx +++ b/src/buildblock/extend_projdata.cxx @@ -30,9 +30,9 @@ START_NAMESPACE_STIR +#if STIR_VERSION < 060000 namespace detail { -#if STIR_VERSION < 060000 /* This function takes symmetries in the sinogram space into account to find data in the negative segment if necessary. However, it needs testing if it would work for non-direct sinograms. diff --git a/src/buildblock/find_STIR_config.cxx b/src/buildblock/find_STIR_config.cxx index 77b290a08..f1bbf9c8a 100644 --- a/src/buildblock/find_STIR_config.cxx +++ b/src/buildblock/find_STIR_config.cxx @@ -26,12 +26,12 @@ START_NAMESPACE_STIR std::string find_STIR_config_file(const std::string& filename){ std::string dir; - dir = get_STIR_config_dir(); //STIR_CONFIG_DIR; + dir = get_STIR_config_dir(); // TODO this might be dangerous on Windows but seems to work const std::string name = (dir+"/"+filename); std::ifstream file(name); if (!file) - error("find_STIR_config_file could not open "+name); + error("find_STIR_config_file could not open " + name + "\nYou can set the STIR_CONFIG_DIR environment variable to help."); if (file.peek() == std::ifstream::traits_type::eof()) error("find_STIR_config_file error opening file for reading (non-existent or empty file). Filename:\n'" + name + "'"); return name; diff --git a/src/data_buildblock/CMakeLists.txt b/src/data_buildblock/CMakeLists.txt index de99dbf5a..e04f5c452 100644 --- a/src/data_buildblock/CMakeLists.txt +++ b/src/data_buildblock/CMakeLists.txt @@ -35,4 +35,4 @@ endif() include(stir_lib_target) -target_link_libraries(${dir} buildblock) +target_link_libraries(${dir} PUBLIC buildblock) diff --git a/src/display/CMakeLists.txt b/src/display/CMakeLists.txt index dc7140814..05d1a5147 100644 --- a/src/display/CMakeLists.txt +++ b/src/display/CMakeLists.txt @@ -46,7 +46,7 @@ include(stir_lib_target) if( "${GRAPHICS}" STREQUAL "X") find_package(Curses REQUIRED) INCLUDE_DIRECTORIES(${X11_INCLUDE_DIR} ${CURSES_INCLUDE_DIR}) - target_link_libraries(${dir} ${X11_LIBRARIES} ${CURSES_LIBRARY}) + target_link_libraries(${dir} PUBLIC ${X11_LIBRARIES} ${CURSES_LIBRARY}) endif() -target_link_libraries(${dir} buildblock) +target_link_libraries(${dir} PUBLIC buildblock) diff --git a/src/eval_buildblock/CMakeLists.txt b/src/eval_buildblock/CMakeLists.txt index 15dd6ac41..2f640cf30 100644 --- a/src/eval_buildblock/CMakeLists.txt +++ b/src/eval_buildblock/CMakeLists.txt @@ -15,6 +15,6 @@ set(${dir_LIB_SOURCES} include(stir_lib_target) -target_link_libraries(eval_buildblock buildblock ) +target_link_libraries(eval_buildblock PUBLIC buildblock ) diff --git a/src/experimental/buildblock/CMakeLists.txt b/src/experimental/buildblock/CMakeLists.txt index d68f8aa9b..519a6461e 100644 --- a/src/experimental/buildblock/CMakeLists.txt +++ b/src/experimental/buildblock/CMakeLists.txt @@ -8,5 +8,5 @@ set(${dir_LIB_SOURCES} include(stir_lib_target) -target_link_libraries(${dir} buildblock) +target_link_libraries(${dir} PUBLIC buildblock) diff --git a/src/experimental/listmode/CMakeLists.txt b/src/experimental/listmode/CMakeLists.txt index 30ae63c8e..845f467f7 100644 --- a/src/experimental/listmode/CMakeLists.txt +++ b/src/experimental/listmode/CMakeLists.txt @@ -16,4 +16,4 @@ endif() include(stir_lib_target) -target_link_libraries(local_listmode_buildblock listmode_buildblock ) +target_link_libraries(local_listmode_buildblock PUBLIC listmode_buildblock ) diff --git a/src/experimental/motion/CMakeLists.txt b/src/experimental/motion/CMakeLists.txt index c9bc2989e..21e283408 100644 --- a/src/experimental/motion/CMakeLists.txt +++ b/src/experimental/motion/CMakeLists.txt @@ -11,5 +11,5 @@ set(${dir_LIB_SOURCES} include(stir_lib_target) -target_link_libraries(${dir} numerics_buildblock local_buildblock ) +target_link_libraries(${dir} PUBLIC numerics_buildblock local_buildblock ) diff --git a/src/experimental/motion_utilities/CMakeLists.txt b/src/experimental/motion_utilities/CMakeLists.txt index 7397b263e..58299a07f 100644 --- a/src/experimental/motion_utilities/CMakeLists.txt +++ b/src/experimental/motion_utilities/CMakeLists.txt @@ -22,7 +22,7 @@ set(${dir_EXE_SOURCES} #include(stir_exe_targets) foreach(executable ${${dir_EXE_SOURCES}}) add_executable(${executable} ${executable} ${STIR_IO_REGISTRIES} ) - target_link_libraries(${executable} buildblock IO buildblock local_motion_buildblock buildblock IO buildblock listmode_buildblock display) + target_link_libraries(${executable} PUBLIC buildblock IO buildblock local_motion_buildblock buildblock IO buildblock listmode_buildblock display) SET_PROPERTY(TARGET ${executable} PROPERTY FOLDER "Executables") endforeach() diff --git a/src/experimental/recon_buildblock/CMakeLists.txt b/src/experimental/recon_buildblock/CMakeLists.txt index 668004ebc..7abfbca26 100644 --- a/src/experimental/recon_buildblock/CMakeLists.txt +++ b/src/experimental/recon_buildblock/CMakeLists.txt @@ -22,7 +22,7 @@ set(${dir_LIB_SOURCES} include(stir_lib_target) -target_link_libraries(local_recon_buildblock display buildblock recon_buildblock) +target_link_libraries(local_recon_buildblock PUBLIC display buildblock recon_buildblock) diff --git a/src/experimental/utilities/list_TAC_ROI_values.cxx b/src/experimental/utilities/list_TAC_ROI_values.cxx index de1d8a80a..cd739302b 100644 --- a/src/experimental/utilities/list_TAC_ROI_values.cxx +++ b/src/experimental/utilities/list_TAC_ROI_values.cxx @@ -23,7 +23,7 @@ \param output_filename a text file sorted in list form of: | Frame_num | Start Time | End Time | Mean | StdDev | CV | \param data_filename should be in ECAT7 format. \param --CV is used to output the Coefficient of Variation as well. - \Note When ROI_filename.par is not given, the user will be asked for the parameters. + Note: When ROI_filename.par is not given, the user will be asked for the parameters. The .par file has the following format \verbatim diff --git a/src/include/stir/Bin.h b/src/include/stir/Bin.h index f2dd3a803..f39066f3d 100644 --- a/src/include/stir/Bin.h +++ b/src/include/stir/Bin.h @@ -18,6 +18,7 @@ Copyright (C) 2000 PARAPET partners Copyright (C) 2000- 2009, Hammersmith Imanet Ltd Copyright (C) 2016, University of Hull + Copyright (C) 2023, University College London This file is part of STIR. SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license @@ -28,7 +29,7 @@ #define __stir_Bin_H__ -#include "stir/common.h" +#include "stir/ViewgramIndices.h" START_NAMESPACE_STIR @@ -46,43 +47,37 @@ START_NAMESPACE_STIR changes to all the framework, when one set a float value, it has to be as 'x.f' */ -class Bin +class Bin : public ViewgramIndices { + typedef ViewgramIndices base_type; public: - //! default constructor + //! default constructor (leaves most members uninitialised) inline Bin(); - //! A constructor : constructs a bin with value (defaulting to 0) + //! constructs a bin with value (timing_pos defaulting to 0) inline Bin(int segment_num,int view_num, int axial_pos_num, int tangential_pos_num,float bin_value); + //! constructs a bin with value (timing_pos and value defaulting to 0) inline Bin(int segment_num, int view_num, int axial_pos_num, int tangential_pos_num); inline Bin(int segment_num, int view_num, int axial_pos_num, int tangential_pos_num, int timing_pos_num, float bin_value); + //! constructs a bin with value (value defaulting to 0) inline Bin(int segment_num, int view_num, int axial_pos_num, int tangential_pos_num, int timing_pos_num); //!get axial position number inline int axial_pos_num()const; - //! get segmnet number - inline int segment_num()const; //! get tangential position number inline int tangential_pos_num() const; - //! get view number - inline int view_num() const; - //! get timing position number - inline int timing_pos_num() const; //! get time-frame number (1-based) inline int time_frame_num() const; inline int& axial_pos_num(); - inline int& segment_num(); inline int& tangential_pos_num(); - inline int& view_num(); - inline int& timing_pos_num(); inline int& time_frame_num(); //! get an empty copy @@ -109,11 +104,8 @@ class Bin private : // shared_ptr proj_data_info_ptr; - int segment; - int view; int axial_pos; int tangential_pos; - int timing_pos; float bin_value; int time_frame; diff --git a/src/include/stir/Bin.inl b/src/include/stir/Bin.inl index d6bd083b3..f6f618c05 100644 --- a/src/include/stir/Bin.inl +++ b/src/include/stir/Bin.inl @@ -17,6 +17,7 @@ Copyright (C) 2000 PARAPET partners Copyright (C) 2000- 2009, Hammersmith Imanet Ltd Copyright (C) 2016, University of Hull + Copyright (C) 2023, University College London This file is part of STIR. SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license @@ -26,52 +27,38 @@ START_NAMESPACE_STIR -Bin::Bin():segment(0),view(0), - axial_pos(0),tangential_pos(0),timing_pos(0), bin_value(0.0f),time_frame(1) +Bin::Bin() + : axial_pos(0),tangential_pos(0), bin_value(0.0f),time_frame(1) {} - Bin::Bin(int segment_num,int view_num, int axial_pos_num,int tangential_pos_num, float bin_value) - :segment(segment_num),view(view_num), - axial_pos(axial_pos_num),tangential_pos(tangential_pos_num), timing_pos(0), bin_value(bin_value),time_frame(1) - {} + : ViewgramIndices(view_num, segment_num, /*timing_pos_num*/0), + axial_pos(axial_pos_num),tangential_pos(tangential_pos_num), bin_value(bin_value),time_frame(1) +{} Bin::Bin(int segment_num,int view_num, int axial_pos_num,int tangential_pos_num) - :segment(segment_num),view(view_num), - axial_pos(axial_pos_num),tangential_pos(tangential_pos_num),timing_pos(0), bin_value(0.0f),time_frame(1) - {} + : ViewgramIndices(view_num, segment_num, /*timing_pos_num*/0), + axial_pos(axial_pos_num),tangential_pos(tangential_pos_num), bin_value(0.0f),time_frame(1) +{} Bin::Bin(int segment_num,int view_num, int axial_pos_num,int tangential_pos_num, int timing_pos_num, float bin_value) - :segment(segment_num),view(view_num), - axial_pos(axial_pos_num),tangential_pos(tangential_pos_num), timing_pos(timing_pos_num), bin_value(bin_value),time_frame(1) - {} + : ViewgramIndices(view_num, segment_num, timing_pos_num), + axial_pos(axial_pos_num),tangential_pos(tangential_pos_num), bin_value(bin_value),time_frame(1) +{} Bin::Bin(int segment_num,int view_num, int axial_pos_num,int tangential_pos_num, int timing_pos_num) - :segment(segment_num),view(view_num), - axial_pos(axial_pos_num),tangential_pos(tangential_pos_num), timing_pos(timing_pos_num), bin_value(0.0f), time_frame(1) - {} - + : ViewgramIndices(view_num, segment_num, timing_pos_num), + axial_pos(axial_pos_num),tangential_pos(tangential_pos_num), bin_value(0.0f), time_frame(1) +{} int Bin:: axial_pos_num()const { return axial_pos;} -int - Bin::segment_num()const -{return segment;} - int Bin::tangential_pos_num() const { return tangential_pos;} -int - Bin::view_num() const -{ return view;} - -int -Bin::timing_pos_num() const -{ return timing_pos; } - int Bin:: time_frame_num() const {return time_frame;} @@ -80,24 +67,12 @@ int& Bin::axial_pos_num() { return axial_pos;} -int& - Bin:: segment_num() -{return segment;} - int& Bin::tangential_pos_num() { return tangential_pos;} int& - Bin:: view_num() -{ return view;} - -int& -Bin:: timing_pos_num() -{ return timing_pos;} - -int& -Bin:: time_frame_num() + Bin:: time_frame_num() {return time_frame;} #if 0 @@ -135,9 +110,8 @@ bool Bin::operator==(const Bin& bin2) const { return - segment == bin2.segment && view == bin2.view && + base_type::operator==(bin2) && axial_pos == bin2.axial_pos && tangential_pos == bin2.tangential_pos && - timing_pos == bin2.timing_pos && time_frame == bin2.time_frame && bin_value == bin2.bin_value; } diff --git a/src/include/stir/DataSymmetriesForViewSegmentNumbers.h b/src/include/stir/DataSymmetriesForViewSegmentNumbers.h index a3e6cf0e1..c557b4c9d 100644 --- a/src/include/stir/DataSymmetriesForViewSegmentNumbers.h +++ b/src/include/stir/DataSymmetriesForViewSegmentNumbers.h @@ -22,13 +22,11 @@ #ifndef __DataSymmetriesForViewSegmentNumbers_H__ #define __DataSymmetriesForViewSegmentNumbers_H__ -#include "stir/common.h" +#include "stir/ViewSegmentNumbers.h" #include START_NAMESPACE_STIR -class ViewSegmentNumbers; - #if 0 class ViewSegmentIndexRange; #endif @@ -46,6 +44,8 @@ class ViewSegmentIndexRange; The class mainly defines members to find \c basic ViewSegmentNumbers. These form a 'basis' for all ViewSegmentNumbers in the sense that all ViewSegmentNumbers can be obtained by using symmetry operations on the 'basic' ones. + + \par Warning: This class wil be renamed/revised to work with \c ViewgramIndices instead. */ class DataSymmetriesForViewSegmentNumbers { @@ -74,12 +74,26 @@ class DataSymmetriesForViewSegmentNumbers virtual void get_related_view_segment_numbers(std::vector&, const ViewSegmentNumbers& v_s) const = 0; +#if 0 + // not yet, as would need copying of vector + //! fills in a vector with all the view/segments that are related to 'v_s' (including itself) + virtual std::vector + get_related_view_segment_numbers(const ViewgramIndices& ind) const + { + } +#endif + //! returns the number of view_segment_numbers related to 'v_s' /*! The default implementation is in terms of get_related_view_segment_numbers, which will be slow of course */ virtual int num_related_view_segment_numbers(const ViewSegmentNumbers& v_s) const; + std::size_t num_related_viewgram_indices(const ViewgramIndices& ind) const + { + return static_cast(num_related_view_segment_numbers(ind)); + } + /*! \brief given an arbitrary view/segment, find the basic view/segment sets 'v_s' to the corresponding 'basic' view/segment and returns true if diff --git a/src/include/stir/DetectorCoordinateMap.h b/src/include/stir/DetectorCoordinateMap.h index 359d6be51..ecc79bf19 100644 --- a/src/include/stir/DetectorCoordinateMap.h +++ b/src/include/stir/DetectorCoordinateMap.h @@ -27,6 +27,9 @@ #include #include #include +#ifdef STIR_OPENMP +#include +#endif #include "stir/CartesianCoordinate3D.h" #include "stir/DetectionPosition.h" @@ -44,7 +47,7 @@ class Succeeded; ring,detector,(layer,)x,y,z An empty line will terminate the reading at that line. - Optionally LOR end-points can be randomly displaced using a Gaussian distribution with standard deviation \sigma (in mm). + Optionally LOR end-points can be randomly displaced using a Gaussian distribution with standard deviation sigma (in mm). */ class DetectorCoordinateMap { @@ -65,14 +68,16 @@ class DetectorCoordinateMap //! Constructor calls read_detectormap_from_file( filename ). DetectorCoordinateMap(const std::string& filename, double sigma = 0.0) : - sigma(sigma), - distribution(0.0, sigma) - { read_detectormap_from_file( filename ); } + DetectorCoordinateMap(sigma) + { + read_detectormap_from_file(filename); + } //! Constructor calls set_detector_map(coord_map). DetectorCoordinateMap(const det_pos_to_coord_type& coord_map, double sigma = 0.0) : - sigma(sigma), - distribution(0.0, sigma) - { set_detector_map( coord_map ); } + DetectorCoordinateMap(sigma) + { + set_detector_map( coord_map ); + } //! Reads map from file and stores it. void read_detectormap_from_file( const std::string& filename ); @@ -83,14 +88,23 @@ class DetectorCoordinateMap stir::DetectionPosition<> get_det_pos_for_index(const stir::DetectionPosition<>& index) const { return input_index_to_det_pos.at(index); - } + } //! Returns a cartesian coordinate given a detection position. stir::CartesianCoordinate3D get_coordinate_for_det_pos( const stir::DetectionPosition<>& det_pos ) const { auto coord = det_pos_to_coord.at(det_pos); - coord.x() += static_cast(distribution(generator)); - coord.y() += static_cast(distribution(generator)); - coord.z() += static_cast(distribution(generator)); + if (sigma == 0.0) + return coord; + +#ifdef STIR_OPENMP + auto thread_id = omp_get_thread_num(); +#else + auto thread_id = 0; +#endif + + coord.x() += static_cast(distributions[thread_id](generators[thread_id])); + coord.y() += static_cast(distributions[thread_id](generators[thread_id])); + coord.z() += static_cast(distributions[thread_id](generators[thread_id])); return coord; } //! Returns a cartesian coordinate given an (unsorted) index. @@ -112,10 +126,20 @@ class DetectorCoordinateMap protected: - DetectorCoordinateMap(double sigma = 0.0) : - sigma(sigma), - distribution(0.0, sigma) - {} + explicit DetectorCoordinateMap(double sigma = 0.0) : + sigma(sigma) + { +#ifdef STIR_OPENMP + generators.resize(omp_get_max_threads()); + distributions.resize(omp_get_max_threads()); + for (auto &distribution : distributions) + distribution = std::normal_distribution(0.0, sigma); +#else + generators.resize(1); + distributions.resize(1); + distributions[0] = std::normal_distribution(0.0, sigma); +#endif + } private: unsigned num_tangential_coords; unsigned num_axial_coords; @@ -128,8 +152,8 @@ class DetectorCoordinateMap stir::DetectionPosition<>> detection_position_map_given_cartesian_coord_keys_2_decimal; const double sigma; - mutable std::default_random_engine generator; - mutable std::normal_distribution distribution; + mutable std::vector generators; + mutable std::vector> distributions; static det_pos_to_coord_type read_detectormap_from_file_help( const std::string& crystal_map_name ); diff --git a/src/include/stir/GeometryBlocksOnCylindrical.h b/src/include/stir/GeometryBlocksOnCylindrical.h index a14f5cea5..fba4b1802 100644 --- a/src/include/stir/GeometryBlocksOnCylindrical.h +++ b/src/include/stir/GeometryBlocksOnCylindrical.h @@ -36,12 +36,12 @@ START_NAMESPACE_STIR \ingroup projdata \brief A helper class to build the crystal map based on scanner info. - \This class builds two maps between cartesian coordinates (z, y, x) - \ and the corresponding detection position (tangential_num, axial_num, radial_num) for each crystal. - \The crystal map, then, is used in ProjDataInfoBlocksOnCylindrical, ProjDataInfoBlocksOnCylindricalNoArcCorr, and CListRecordSAFIR + This class builds two maps between cartesian coordinates (z, y, x) + and the corresponding detection position (tangential_num, axial_num, radial_num) for each crystal. + The crystal map, then, is used in ProjDataInfoBlocksOnCylindrical, ProjDataInfoBlocksOnCylindricalNoArcCorr, and CListRecordSAFIR - \The center of first ring is the center of coordinates. - \Distances are from center to center of crystals. + The center of first ring is the center of coordinates. + Distances are from center to center of crystals. */ diff --git a/src/include/stir/HighResWallClockTimer.h b/src/include/stir/HighResWallClockTimer.h index a72346579..eb1ce6a6b 100644 --- a/src/include/stir/HighResWallClockTimer.h +++ b/src/include/stir/HighResWallClockTimer.h @@ -17,7 +17,8 @@ \author Alexey Zverovich \author PARAPET project - +*/ +/* Modification history: diff --git a/src/include/stir/IO/InterfileHeaderSiemens.h b/src/include/stir/IO/InterfileHeaderSiemens.h index 19d41cd17..66a9da645 100644 --- a/src/include/stir/IO/InterfileHeaderSiemens.h +++ b/src/include/stir/IO/InterfileHeaderSiemens.h @@ -1,24 +1,22 @@ /* - Copyright (C) 2002-2007, Hammersmith Imanet Ltd - Copyright (C) 2013, 2016 University College London + Copyright (C) 2020, 2021, 2023 University College London + Copyright (C) 2018 STFC This file is part of STIR. - SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license + SPDX-License-Identifier: Apache-2.0 See STIR/LICENSE.txt for details */ /*! \file \ingroup InterfileIO + \ingroup ECAT \brief This file declares the classes stir::InterfileHeaderSiemens, - stir::InterfileImageHeader, stir::InterfilePDFSHeader + stir::InterfileRawDataHeaderSiemens, stir::InterfilePDFSHeaderSiemens, + stir::InterfileListmodeHeaderSiemens, stir:InterfileNormHeaderSiemens \author Kris Thielemans - \author Sanida Mustafovic - \author PARAPET project - - See http://stir.sourceforge.net for a description of the full - proposal for Interfile headers for 3D PET. + \author Edoardo Pasca */ @@ -39,6 +37,7 @@ class ProjDataInfo; \brief a class for Interfile keywords (and parsing) common to all types of data \ingroup InterfileIO + \ingroup ECAT */ class InterfileHeaderSiemens : public InterfileHeader { @@ -53,6 +52,19 @@ class InterfileHeaderSiemens : public InterfileHeader protected: // Returns false if OK, true if not. virtual bool post_processing(); + //! ignore multiple GMT times + /*! + Siemens uses keywords like + \verbatim + %study date (yyyy:mm:dd") := ... + %study time (hh:mm:ss GMT+00:00) := ... + \endvarbatim + You can ignore this for all (?) time zones by using + \code + ignore_Siemens_date_and_time_keys("%study"); + \endcode + */ + void ignore_Siemens_date_and_time_keys(const std::string& keyword); private: @@ -80,6 +92,7 @@ class InterfileHeaderSiemens : public InterfileHeader /*! \brief a class for Interfile keywords (and parsing) specific to images \ingroup InterfileIO + \ingroup ECAT */ class InterfileImageHeader : public InterfileHeaderSiemens { @@ -101,8 +114,9 @@ class InterfileImageHeader : public InterfileHeaderSiemens /*! \brief a class for Interfile keywords (and parsing) specific to -Siemens PET projection or list mode data +Siemens PET projection, list mode data or norm data \ingroup InterfileIO +\ingroup ECAT */ class InterfileRawDataHeaderSiemens : public InterfileHeaderSiemens { @@ -145,8 +159,9 @@ class InterfileRawDataHeaderSiemens : public InterfileHeaderSiemens /*! \brief a class for Interfile keywords (and parsing) specific to -projection data (i.e. ProjDataFromStream) +projection data (i.e. ProjDataFromStream) for Siemens PET scanners \ingroup InterfileIO +\ingroup ECAT */ class InterfilePDFSHeaderSiemens : public InterfileRawDataHeaderSiemens { @@ -181,8 +196,9 @@ class InterfilePDFSHeaderSiemens : public InterfileRawDataHeaderSiemens /*! \brief a class for Interfile keywords (and parsing) specific to -Siemesn listmode data (in PETLINK format) +Siemens listmode data (in PETLINK format) \ingroup InterfileIO +\ingroup ECAT */ class InterfileListmodeHeaderSiemens : public InterfileRawDataHeaderSiemens { @@ -211,6 +227,34 @@ class InterfileListmodeHeaderSiemens : public InterfileRawDataHeaderSiemens }; +/*! +\brief a class for Interfile keywords (and parsing) specific to +Siemens (component-based) normalisation data +\ingroup InterfileIO +\ingroup ECAT +*/ +class InterfileNormHeaderSiemens : public InterfileRawDataHeaderSiemens + { + public: + InterfileNormHeaderSiemens(); + + protected: + + //! Returns false if OK, true if not. + virtual bool post_processing() override; + + public: + float calib_factor; + float cross_calib_factor; + int num_buckets; + + private: + int num_components; + std::vector> number_of_dimensions; + void read_num_components(); + + }; + END_NAMESPACE_STIR #endif // __stir_InterfileHeaderSiemens_H__ diff --git a/src/include/stir/KeyParser.h b/src/include/stir/KeyParser.h index 1bc022303..175c33765 100644 --- a/src/include/stir/KeyParser.h +++ b/src/include/stir/KeyParser.h @@ -110,7 +110,9 @@ public : \brief A class to parse Interfile headers Currently, Interfile 3.3 parsing rules are hard-coded, with - some extensions. + some extensions. Note that some functions such as `get_keyword()` are `virtual` + such that a derived class could use non-Interfile parsing (but this might need + more work). KeyParser reads input line by line and parses each line separately. It allows for '\\r' at the end of the line (as in files @@ -349,16 +351,14 @@ protected : //! convert 'rough' keyword into a standardised form - /*! \todo Implementation note: this function is non-static such that it can - be overloaded. Probably a template with a function object would be - better. */ + /*! Calls standardise_interfile_keyword(). +. */ virtual std::string standardise_keyword(const std::string& keyword) const; //! gets a keyword from a string - /*! Implementation note: this function is non-static as it uses - standardise_keyword(). + /*! Find `:=` or `[`. Note that the returned keyword is not standardised yet. */ virtual std::string get_keyword(const std::string&) const; diff --git a/src/include/stir/ProjData.h b/src/include/stir/ProjData.h index 28c733414..2cd068b6a 100644 --- a/src/include/stir/ProjData.h +++ b/src/include/stir/ProjData.h @@ -2,7 +2,7 @@ Copyright (C) 2000 PARAPET partners Copyright (C) 2000- 2012, Hammersmith Imanet Ltd Copyright (C) 2016, 2017, University of Hull - Copyright (C) 2013, 2015-2017, 2020, University College London + Copyright (C) 2013, 2015-2017, 2020, 2023 University College London This file is part of STIR. SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license @@ -30,6 +30,9 @@ #include "stir/Succeeded.h" #include "stir/SegmentBySinogram.h" #include "stir/SegmentByView.h" +#include "stir/SegmentIndices.h" +#include "stir/ViewgramIndices.h" +#include "stir/SinogramIndices.h" //#include @@ -44,7 +47,6 @@ template class SegmentBySinogram; template class SegmentByView; template class Viewgram; template class Sinogram; -class ViewSegmentNumbers; class Succeeded; class ProjDataInMemory; //class ExamInfo; @@ -121,18 +123,32 @@ class ProjData : public ExamData inline shared_ptr get_proj_data_info_sptr() const; //! Get viewgram + /*! + \deprecated Use get_viewgram(const ViewgramIndices&) instead. + */ virtual Viewgram get_viewgram(const int view, const int segment_num, const bool make_num_tangential_poss_odd = false, const int timing_pos = 0) const = 0; + //! Get viewgram + inline Viewgram + get_viewgram(const ViewgramIndices&) const; + //! Set viewgram virtual Succeeded set_viewgram(const Viewgram&) = 0; //! Get sinogram + /*! + \deprecated Use get_sinogram(const SinogramIndices&) instead . + */ virtual Sinogram get_sinogram(const int ax_pos_num, const int segment_num, const bool make_num_tangential_poss_odd = false, const int timing_pos = 0) const = 0; + //! Get sinogram + inline Sinogram + get_sinogram(const SinogramIndices&) const; + //! Set sinogram virtual Succeeded set_sinogram(const Sinogram&) = 0; @@ -144,32 +160,72 @@ class ProjData : public ExamData get_subset(const std::vector& views) const; //! Get empty viewgram + Viewgram get_empty_viewgram(const ViewgramIndices&) const; + + //! Get empty viewgram + /*! + \deprecated Use get_viewgram(const ViewgramIndices&) instead. + */ Viewgram get_empty_viewgram(const int view, const int segment_num, const bool make_num_tangential_poss_odd = false, const int timing_pos = 0) const; //! Get empty_sinogram + Sinogram + get_empty_sinogram(const SinogramIndices&) const; + + //! Get empty_sinogram + /*! + \deprecated Use get_sinogram(const SinogramIndices&) instead . + */ Sinogram get_empty_sinogram(const int ax_pos_num, const int segment_num, const bool make_num_tangential_poss_odd = false, const int timing_pos = 0) const; - //! Get empty segment sino - SegmentByView + //! Get empty segment by view + SegmentByView + get_empty_segment_by_view(const SegmentIndices&) const; + //! Get empty segment by sino + SegmentBySinogram + get_empty_segment_by_sinogram(const SegmentIndices&) const; + //! Get empty segment view + /*! + \deprecated Use get_empty_segment_by_sinogram(const SegmentIndices&) instead . + */ + SegmentByView get_empty_segment_by_view(const int segment_num, const bool make_num_tangential_poss_odd = false, const int timing_pos = 0) const; - //! Get empty segment view - SegmentBySinogram + //! Get empty segment sino + /*! + \deprecated Use get_empty_segment_by_sinogram(const SegmentIndices&) instead . + */ + SegmentBySinogram get_empty_segment_by_sinogram(const int segment_num, const bool make_num_tangential_poss_odd = false, const int timing_pos = 0) const; - //! Get segment by sinogram + /*! + \deprecated Use get_segment_by_sinogram(const SegmentIndices&) instead. + */ virtual SegmentBySinogram get_segment_by_sinogram(const int segment_num, const int timing_pos = 0) const; + + //! Get segment by sinogram + inline SegmentBySinogram + get_segment_by_sinogram(const SegmentIndices&) const; + //! Get segment by view - virtual SegmentByView + /*! + \deprecated Use get_segment_by_view(const SegmentIndices&) instead. + */ + virtual SegmentByView get_segment_by_view(const int segment_num, const int timing_pos = 0) const; + + //! Get segment by view + inline SegmentByView + get_segment_by_view(const SegmentIndices&) const; + //! Set segment by sinogram virtual Succeeded set_segment(const SegmentBySinogram&); @@ -178,8 +234,9 @@ class ProjData : public ExamData set_segment(const SegmentByView&); //! Get related viewgrams + // TODOTOF remove timing_pos arg virtual RelatedViewgrams - get_related_viewgrams(const ViewSegmentNumbers&, + get_related_viewgrams(const ViewgramIndices&, const shared_ptr&, const bool make_num_tangential_poss_odd = false, const int timing_pos = 0) const; @@ -191,9 +248,8 @@ class ProjData : public ExamData //! Get empty related viewgrams, where the symmetries_ptr specifies the symmetries to use RelatedViewgrams - get_empty_related_viewgrams(const ViewSegmentNumbers& view_segmnet_num, - //const int view_num, const int segment_num, - const shared_ptr& symmetries_ptr, + get_empty_related_viewgrams(const ViewgramIndices& viewgram_indices, + const shared_ptr& symmetries_ptr, const bool make_num_tangential_poss_odd = false, const int timing_pos = 0) const; diff --git a/src/include/stir/ProjData.inl b/src/include/stir/ProjData.inl index 6f028540b..59cc165c9 100644 --- a/src/include/stir/ProjData.inl +++ b/src/include/stir/ProjData.inl @@ -11,7 +11,7 @@ /* Copyright (C) 2000 PARAPET partners Copyright (C) 2000-2009, Hammersmith Imanet Ltd - Copyright (C) 2013, 2015 University College London + Copyright (C) 2013, 2015, 2023 University College London Copyright (C) 2016, University of Hull This file is part of STIR. @@ -24,6 +24,30 @@ START_NAMESPACE_STIR +SegmentBySinogram +ProjData::get_segment_by_sinogram(const SegmentIndices& si) const +{ + return this->get_segment_by_sinogram(si.segment_num()); +} + +SegmentByView +ProjData::get_segment_by_view(const SegmentIndices& si) const +{ + return this->get_segment_by_view(si.segment_num()); +} + +Viewgram +ProjData::get_viewgram(const ViewgramIndices& vi) const +{ + return this->get_viewgram(vi.view_num(), vi.segment_num(), false, vi.timing_pos_num()); +} + +Sinogram +ProjData::get_sinogram(const SinogramIndices& vi) const +{ + return this->get_sinogram(vi.axial_pos_num(), vi.segment_num(), false, vi.timing_pos_num()); +} + shared_ptr ProjData::get_proj_data_info_sptr() const { diff --git a/src/include/stir/ProjDataInfo.h b/src/include/stir/ProjDataInfo.h index 8b0b27953..d7d232b7f 100644 --- a/src/include/stir/ProjDataInfo.h +++ b/src/include/stir/ProjDataInfo.h @@ -5,7 +5,7 @@ Copyright (C) 2000 - 2011-10-14, Hammersmith Imanet Ltd Copyright (C) 2011-07-01 - 2011, Kris Thielemans Copyright (C) 2016-17, University of Hull - Copyright (C) 2017-2018, 2020, 2022, University College London + Copyright (C) 2017-2018, 2020, 2022, 2023 University College London This file is part of STIR. SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license @@ -27,6 +27,9 @@ #ifndef __stir_ProjDataInfo_H__ #define __stir_ProjDataInfo_H__ +#include "stir/SegmentIndices.h" +#include "stir/ViewgramIndices.h" +#include "stir/SinogramIndices.h" #include "stir/VectorWithOffset.h" #include "stir/Scanner.h" #include "stir/shared_ptr.h" @@ -405,25 +408,40 @@ class ProjDataInfo //@{ //! Get empty viewgram + Viewgram get_empty_viewgram(const ViewgramIndices&) const; + //! Get empty_sinogram + Sinogram get_empty_sinogram(const SinogramIndices&) const; + //! Get empty segment sino + SegmentByView get_empty_segment_by_view(const SegmentIndices&) const; + //! Get empty segment view + SegmentBySinogram get_empty_segment_by_sinogram(const SegmentIndices&) const; + + //! Get empty viewgram + /*! \deprecated */ Viewgram get_empty_viewgram(const int view_num, const int segment_num, const bool make_num_tangential_poss_odd = false, const int timing_pos_num = 0) const; //! Get empty_sinogram + /*! \deprecated */ Sinogram get_empty_sinogram(const int ax_pos_num, const int segment_num, const bool make_num_tangential_poss_odd = false, const int timing_pos_num = 0) const; //! Get empty segment sino + /*! \deprecated */ SegmentByView get_empty_segment_by_view(const int segment_num, const bool make_num_tangential_poss_odd = false, const int timing_pos_num = 0) const; //! Get empty segment view + /*! \deprecated */ SegmentBySinogram get_empty_segment_by_sinogram(const int segment_num, const bool make_num_tangential_poss_odd = false, const int timing_pos_num = 0) const; //! Get empty related viewgrams, where the symmetries_ptr specifies the symmetries to use - RelatedViewgrams get_empty_related_viewgrams(const ViewSegmentNumbers&, + /*! make_num_tangential_poss_odd has to be \c false */ + //TODOTOF + RelatedViewgrams get_empty_related_viewgrams(const ViewgramIndices&, const shared_ptr&, const bool make_num_tangential_poss_odd = false, const int timing_pos_num = 0) const; diff --git a/src/include/stir/RelatedViewgrams.h b/src/include/stir/RelatedViewgrams.h index 0649cd88b..f6c632db9 100644 --- a/src/include/stir/RelatedViewgrams.h +++ b/src/include/stir/RelatedViewgrams.h @@ -88,8 +88,11 @@ class RelatedViewgrams /*! see DataSymmetriesForViewSegmentNumbers for definition of 'basic' */ inline int get_basic_timing_pos_num() const; //! get 'basic' view_segment_num + /*! \deprecated Use get_basic_viewgram_indices() instead. */ + inline ViewgramIndices get_basic_view_segment_num() const; + //! get 'basic' viewgram indices /*! see DataSymmetriesForViewSegmentNumbers for definition of 'basic' */ - inline ViewSegmentNumbers get_basic_view_segment_num() const; + inline ViewgramIndices get_basic_viewgram_indices() const; //! returns the number of viewgrams in this object inline int get_num_viewgrams() const; diff --git a/src/include/stir/RelatedViewgrams.inl b/src/include/stir/RelatedViewgrams.inl index abac77e9f..b3adeeb8f 100644 --- a/src/include/stir/RelatedViewgrams.inl +++ b/src/include/stir/RelatedViewgrams.inl @@ -79,10 +79,21 @@ int RelatedViewgrams::get_basic_timing_pos_num() const } template -ViewSegmentNumbers RelatedViewgrams:: +ViewgramIndices +RelatedViewgrams:: +get_basic_viewgram_indices() const +{ + assert(viewgrams.size()>0); + check_state(); + return viewgrams[0].get_viewgram_indices(); +} + +template +ViewgramIndices +RelatedViewgrams:: get_basic_view_segment_num() const { - return ViewSegmentNumbers(get_basic_view_num(), get_basic_segment_num()); + return this->get_basic_viewgram_indices(); } template diff --git a/src/include/stir/Segment.h b/src/include/stir/Segment.h index 17005b195..72e5678ad 100644 --- a/src/include/stir/Segment.h +++ b/src/include/stir/Segment.h @@ -1,6 +1,7 @@ /* Copyright (C) 2000 PARAPET partners Copyright (C) 2000-2012 Hammersmith Imanet Ltd + Copyright (C) 2023, University College London This file is part of STIR. SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license @@ -21,6 +22,7 @@ #include "stir/ProjDataInfo.h" +#include "stir/SegmentIndices.h" #include "stir/shared_ptr.h" START_NAMESPACE_STIR @@ -33,7 +35,7 @@ template class Viewgram; \ingroup projdata This stores a subset of the data accessible via a ProjData object, - where the segment_num is fixed. + where the SegmentIndices are fixed. At the moment, 2 'storage modes' are supported (and implemented as derived classes). @@ -60,6 +62,7 @@ class Segment get_proj_data_info_sptr() const; virtual StorageOrder get_storage_order() const = 0; + inline SegmentIndices get_segment_indices() const; //! Get the segment number inline int get_segment_num() const; //! Get the timing position index @@ -116,11 +119,9 @@ class Segment protected: shared_ptr proj_data_info_sptr; - int segment_num; - int timing_pos_num; + SegmentIndices _indices; - inline Segment(const shared_ptr& proj_data_info_ptr_v,const int s_num, const int t_num = 0); - + inline Segment(const shared_ptr& proj_data_info_sptr_v,const SegmentIndices&); }; END_NAMESPACE_STIR diff --git a/src/include/stir/Segment.inl b/src/include/stir/Segment.inl index 2cad99189..6623fe622 100644 --- a/src/include/stir/Segment.inl +++ b/src/include/stir/Segment.inl @@ -3,6 +3,7 @@ /* Copyright (C) 2000 PARAPET partners Copyright (C) 2000- 2007, IRSL + Copyright (C) 2023, University College London This file is part of STIR. SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license @@ -25,22 +26,26 @@ START_NAMESPACE_STIR template Segment:: -Segment( const shared_ptr& proj_data_info_ptr_v,const int s_num, const int t_num) +Segment( const shared_ptr& proj_data_info_sptr_v,const SegmentIndices& ind) : - proj_data_info_sptr(proj_data_info_ptr_v), - segment_num(s_num), - timing_pos_num(t_num) + proj_data_info_sptr(proj_data_info_sptr_v), + _indices(ind) {} +template +SegmentIndices +Segment:: get_segment_indices() const +{ return _indices; } + template int Segment:: get_segment_num() const -{ return segment_num; } +{ return _indices.segment_num(); } template int Segment:: get_timing_pos_num() const -{ return timing_pos_num; } +{ return _indices.timing_pos_num(); } template shared_ptr diff --git a/src/include/stir/SegmentBySinogram.h b/src/include/stir/SegmentBySinogram.h index c44190ade..a3a72143b 100644 --- a/src/include/stir/SegmentBySinogram.h +++ b/src/include/stir/SegmentBySinogram.h @@ -3,6 +3,7 @@ /* Copyright (C) 2000 PARAPET partners Copyright (C) 2000- 2012, Hammersmith Imanet Ltd + Copyright (C) 2023, University College London This file is part of STIR. SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license @@ -38,7 +39,7 @@ template class SegmentByView; /*! \ingroup projdata - \brief A class for storing (3d) projection data with a fixed segment_num. + \brief A class for storing (3d) projection data with fixed SegmentIndices. Storage order is as follows: \code @@ -56,11 +57,26 @@ class SegmentBySinogram : public Segment, public Array<3,elemT> typedef typename Segment::StorageOrder StorageOrder; //! Constructor that sets the data to a given 3d Array + SegmentBySinogram(const Array<3,elemT>& v, + const shared_ptr& proj_data_info_ptr_v, + const SegmentIndices& ind); + + //! Constructor that sets sizes via the ProjDataInfo object, initialising data to 0 + SegmentBySinogram(const shared_ptr& proj_data_info_ptr_v, + const SegmentIndices& ind); + + //! Constructor that sets the data to a given 3d Array + /*! + \deprecated Use version with SegmentIndices instead + */ SegmentBySinogram(const Array<3,elemT>& v, const shared_ptr& proj_data_info_ptr_v, const int segment_num, const int timing_pos_num = 0); //! Constructor that sets sizes via the ProjDataInfo object, initialising data to 0 + /*! + \deprecated Use version with SegmentIndices instead + */ SegmentBySinogram(const shared_ptr& proj_data_info_ptr_v, const int segment_num, const int timing_pos_num = 0); diff --git a/src/include/stir/SegmentBySinogram.inl b/src/include/stir/SegmentBySinogram.inl index 6118ddb3a..3b6e128a2 100644 --- a/src/include/stir/SegmentBySinogram.inl +++ b/src/include/stir/SegmentBySinogram.inl @@ -3,6 +3,7 @@ /* Copyright (C) 2000 PARAPET partners Copyright (C) 2000- 2007, Hammersmith Imanet Ltd + Copyright (C) 2023, University College London This file is part of STIR. SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license diff --git a/src/include/stir/SegmentByView.h b/src/include/stir/SegmentByView.h index 4358ac37c..f4381c5ec 100644 --- a/src/include/stir/SegmentByView.h +++ b/src/include/stir/SegmentByView.h @@ -3,6 +3,7 @@ /* Copyright (C) 2000 PARAPET partners Copyright (C) 2000- 2007, Hammersmith Imanet Ltd + Copyright (C) 2023, University College London This file is part of STIR. SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license @@ -37,7 +38,7 @@ template class Sinogram; /*! \ingroup projdata - \brief A class for storing (3d) projection data with a fixed segment_num. + \brief A class for storing (3d) projection data with fixed SegmentIndices. Storage order is as follows: \code @@ -55,16 +56,30 @@ template class SegmentByView : public Segment, public Ar typedef typename Segment::StorageOrder StorageOrder; //! Constructor that sets the data to a given 3d Array + SegmentByView(const Array<3,elemT>& v, + const shared_ptr& proj_data_info_sptr, + const SegmentIndices&); + + //! Constructor that sets sizes via the ProjDataInfo object, initialising data to 0 + SegmentByView(const shared_ptr& proj_data_info_sptr, + const SegmentIndices&); + + //! Constructor that sets the data to a given 3d Array + /*! + \deprecated Use version with SegmentIndices instead + */ SegmentByView(const Array<3,elemT>& v, const shared_ptr& proj_data_info_ptr, const int segment_num, const int timing_pos_num = 0); //! Constructor that sets sizes via the ProjDataInfo object, initialising data to 0 + /*! + \deprecated Use version with SegmentIndices instead + */ SegmentByView(const shared_ptr& proj_data_info_ptr, const int segment_num, const int timing_pos_num = 0); - //! Conversion from 1 storage order to the other SegmentByView(const SegmentBySinogram& ); diff --git a/src/include/stir/SegmentByView.inl b/src/include/stir/SegmentByView.inl index 0f89b7904..80341d498 100644 --- a/src/include/stir/SegmentByView.inl +++ b/src/include/stir/SegmentByView.inl @@ -16,6 +16,7 @@ /* Copyright (C) 2000 PARAPET partners Copyright (C) 2000- 2011, Hammersmtih Imanet Ltd + Copyright (C) 2023, University College London This file is part of STIR. SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license diff --git a/src/include/stir/SegmentIndices.h b/src/include/stir/SegmentIndices.h new file mode 100644 index 000000000..ca433f1f5 --- /dev/null +++ b/src/include/stir/SegmentIndices.h @@ -0,0 +1,70 @@ +// +// +/*! + \file + \ingroup projdata + + \brief Definition of class stir::SegmentIndices + + \author Kris Thielemans + +*/ +/* + Copyright (C) 2023, University College London + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +#ifndef __stir_SegmentIndices_h__ +#define __stir_SegmentIndices_h__ + +#include "stir/common.h" + +START_NAMESPACE_STIR + +/*! + \brief A very simple class to store segment numbers and any other + indices that define a segment + \ingroup projdata +*/ +class SegmentIndices +{ +public: + //! constructor segment number as arguments + explicit inline SegmentIndices(const int segment_num = 0, const int timing_pos_num = 0); + + //! get segment number for const objects + inline int segment_num() const; + + //! get reference to segment number + inline int& segment_num(); + + //! get TOF index for const objects + inline int timing_pos_num() const; + + //! get reference to TOF index + inline int& timing_pos_num(); + + //! comparison operator, only useful for sorting + /*! In future, there will be multiple indices, and order will then be based as in + (0,1) < (0,-1) < (1,1) ... + */ + inline bool operator<(const SegmentIndices& other) const; + + //! test for equality + inline bool operator==(const SegmentIndices& other) const; + inline bool operator!=(const SegmentIndices& other) const; + +private: + int _segment; + int _timing_pos; +}; + +END_NAMESPACE_STIR + +#include "stir/SegmentIndices.inl" + +#endif diff --git a/src/include/stir/SegmentIndices.inl b/src/include/stir/SegmentIndices.inl new file mode 100644 index 000000000..91092c0de --- /dev/null +++ b/src/include/stir/SegmentIndices.inl @@ -0,0 +1,71 @@ +/*! + \file + \ingroup projdata + + \brief inline implementations for class stir::SegmentIndices + + \author Kris Thielemans + \author Sanida Mustafovic + \author PARAPET project + + +*/ +/* + Copyright (C) 2000 PARAPET partners + Copyright (C) 2000- 2009, Hammersmith Imanet Ltd + Copyright (C) 2023, University College London + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license + + See STIR/LICENSE.txt for details +*/ + +START_NAMESPACE_STIR + +SegmentIndices::SegmentIndices(const int segment_num, const int timing_pos_num) + : _segment(segment_num), _timing_pos(timing_pos_num) +{} + +int +SegmentIndices::segment_num() const +{ + return _segment; +} + +int& +SegmentIndices::segment_num() +{ + return _segment; +} + +int +SegmentIndices::timing_pos_num() const +{ + return _timing_pos; +} + +int& +SegmentIndices::timing_pos_num() +{ + return _timing_pos; +} + +bool +SegmentIndices::operator<(const SegmentIndices& other) const +{ + return (_segment < other._segment) || ((_segment == other._segment) && (_timing_pos < other._timing_pos)); +} + +bool +SegmentIndices::operator==(const SegmentIndices& other) const +{ + return (_segment == other._segment) && (_timing_pos == other._timing_pos); +} + +bool +SegmentIndices::operator!=(const SegmentIndices& other) const +{ + return !(*this == other); +} +END_NAMESPACE_STIR diff --git a/src/include/stir/Sinogram.h b/src/include/stir/Sinogram.h index c24533f20..728f6ee87 100644 --- a/src/include/stir/Sinogram.h +++ b/src/include/stir/Sinogram.h @@ -4,6 +4,7 @@ Copyright (C) 2000 PARAPET partners Copyright (C) 2000 - 2007-10-08, Hammersmith Imanet Ltd Copyright (C) 2011-07-01 - 2012, Kris Thielemans + Copyright (C) 2023, University College London This file is part of STIR. SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license @@ -29,6 +30,7 @@ #include "stir/Array.h" #include "stir/ProjDataInfo.h" +#include "stir/SinogramIndices.h" #include "stir/shared_ptr.h" @@ -40,7 +42,7 @@ START_NAMESPACE_STIR \ingroup projdata \brief A class for 2d projection data. - This represents a subset of the full projection. segment_num and axial_pos_num + This represents a subset of the full projection. SegmentIndices and axial_pos_num are fixed. */ @@ -60,14 +62,30 @@ class Sinogram : public Array<2,elemT> #endif public: - //! Construct sinogram from proj_data_info pointer, ring and segment number. Data are set to 0. + //! Construct sinogram from proj_data_info pointe and indices. Data are set to 0. + inline Sinogram(const shared_ptr& proj_data_info_sptr, + const SinogramIndices&); + + //! Construct sinogram with data set to the array. + inline Sinogram(const Array<2,elemT>& p,const shared_ptr& proj_data_info_sptr, + const SinogramIndices&); + + //! Construct sinogram from proj_data_info pointer, axial position and segment number. Data are set to 0. + /*! + \deprecated Use version with SinogramIndices instead. + */ inline Sinogram(const shared_ptr& proj_data_info_ptr, const int ax_pos_num, const int segment_num, const int timing_pos_num = 0); //! Construct sinogram with data set to the array. + /*! + \deprecated Use version with SinogramIndices instead. + */ inline Sinogram(const Array<2,elemT>& p,const shared_ptr& proj_data_info_ptr, const int ax_pos_num, const int segment_num, const int timing_pos_num = 0); - + + //! Get indices + inline SinogramIndices get_sinogram_indices() const; //! Get segment number inline int get_segment_num() const; //! Get number of axial positions @@ -131,9 +149,7 @@ class Sinogram : public Array<2,elemT> private: shared_ptr proj_data_info_ptr; - int axial_pos_num; - int segment_num; - int timing_pos_num; + SinogramIndices _indices; }; diff --git a/src/include/stir/Sinogram.inl b/src/include/stir/Sinogram.inl index 8d2446fa8..d8feda9ce 100644 --- a/src/include/stir/Sinogram.inl +++ b/src/include/stir/Sinogram.inl @@ -3,6 +3,7 @@ /* Copyright (C) 2000 PARAPET partners Copyright (C) 2000- 2007,Hammersmith Imanet Ltd + Copyright (C) 2023, University College London This file is part of STIR. SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license @@ -28,20 +29,27 @@ START_NAMESPACE_STIR +template +SinogramIndices +Sinogram::get_sinogram_indices() const +{ + return this->_indices; +} + template int Sinogram::get_segment_num() const -{ return segment_num; } +{ return this->_indices.segment_num(); } template int Sinogram::get_axial_pos_num() const -{ return axial_pos_num; } +{ return this->_indices.axial_pos_num(); } template int Sinogram::get_timing_pos_num() const -{ return timing_pos_num; } +{ return this->_indices.timing_pos_num(); } template int @@ -79,8 +87,8 @@ template Sinogram Sinogram::get_empty_copy(void) const { - Sinogram copy(proj_data_info_ptr, get_axial_pos_num(), get_segment_num(), get_timing_pos_num()); - return copy; + Sinogram copy(proj_data_info_ptr, get_sinogram_indices()); + return copy; } template @@ -94,16 +102,14 @@ template Sinogram:: Sinogram(const Array<2,elemT>& p, const shared_ptr& pdi_ptr, - const int ax_pos_num, const int s_num, const int t_num) + const SinogramIndices& ind) : - Array<2,elemT>(p), + Array<2,elemT>(p), proj_data_info_ptr(pdi_ptr), - axial_pos_num(ax_pos_num), - segment_num(s_num), - timing_pos_num(t_num) + _indices(ind) { - assert(axial_pos_num <= proj_data_info_ptr->get_max_axial_pos_num(segment_num)); - assert(axial_pos_num >= proj_data_info_ptr->get_min_axial_pos_num(segment_num)); + assert(ind.axial_pos_num() <= proj_data_info_ptr->get_max_axial_pos_num(ind.segment_num())); + assert(ind.axial_pos_num() >= proj_data_info_ptr->get_min_axial_pos_num(ind.segment_num())); // segment_num is already checked by doing get_max_axial_pos_num(s_num) assert( get_min_view_num() == pdi_ptr->get_min_view_num()); @@ -117,21 +123,33 @@ Sinogram(const Array<2,elemT>& p, template Sinogram:: Sinogram(const shared_ptr& pdi_ptr, - const int ax_pos_num, const int s_num, const int t_num) + const SinogramIndices& ind) : Array<2,elemT>(IndexRange2D (pdi_ptr->get_min_view_num(), pdi_ptr->get_max_view_num(), pdi_ptr->get_min_tangential_pos_num(), pdi_ptr->get_max_tangential_pos_num())), proj_data_info_ptr(pdi_ptr), - axial_pos_num(ax_pos_num), - segment_num(s_num), - timing_pos_num(t_num) + _indices(ind) { - assert(axial_pos_num <= proj_data_info_ptr->get_max_axial_pos_num(segment_num)); - assert(axial_pos_num >= proj_data_info_ptr->get_min_axial_pos_num(segment_num)); + assert(ind.axial_pos_num() <= proj_data_info_ptr->get_max_axial_pos_num(ind.segment_num())); + assert(ind.axial_pos_num() >= proj_data_info_ptr->get_min_axial_pos_num(ind.segment_num())); // segment_num is already checked by doing get_max_axial_pos_num(s_num) } +template +Sinogram:: +Sinogram(const Array<2,elemT>& p, + const shared_ptr& pdi_sptr, + const int ax_pos_num, const int s_num, const int t_num) + : Sinogram(p, pdi_sptr, SinogramIndices(ax_pos_num, s_num, t_num)) +{} + +template +Sinogram:: +Sinogram(const shared_ptr& pdi_sptr, + const int ax_pos_num, const int s_num, const int t_num) + : Sinogram(pdi_sptr, SinogramIndices(ax_pos_num, s_num, t_num)) +{} END_NAMESPACE_STIR diff --git a/src/include/stir/SinogramIndices.h b/src/include/stir/SinogramIndices.h new file mode 100644 index 000000000..d7ac00aa0 --- /dev/null +++ b/src/include/stir/SinogramIndices.h @@ -0,0 +1,68 @@ +// +// +/*! + \file + \ingroup projdata + + \brief Definition of class stir::SinogramIndices + + \author Kris Thielemans + +*/ +/* + Copyright (C) 2023, University College London + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + + + +#ifndef __stir_SinogramIndices_h__ +#define __stir_SinogramIndices_h__ + +#include "stir/SegmentIndices.h" + +START_NAMESPACE_STIR + +/*! + \brief A very simple class to store all dincies to get a (2D) Sinogram + \ingroup projdata +*/ +class SinogramIndices : public SegmentIndices +{ + typedef SegmentIndices base_type; +public: + + //! an empty constructor (sets everything to 0) + inline SinogramIndices(); + //! constructor specifying indices + inline SinogramIndices( const int axial_pos_num,const int segment_num, const int timing_pos_num); + + //! get view number for const objects + inline int axial_pos_num() const; + + //! get reference to view number + inline int& axial_pos_num(); + + + //! comparison operator, only useful for sorting + /*! order : (0,1) < (0,-1) < (1,1) ...*/ + inline bool operator<(const SinogramIndices& other) const; + + //! test for equality + inline bool operator==(const SinogramIndices& other) const; + inline bool operator!=(const SinogramIndices& other) const; + +private: + int _axial_pos; + +}; + +END_NAMESPACE_STIR + +#include "stir/SinogramIndices.inl" + +#endif diff --git a/src/include/stir/SinogramIndices.inl b/src/include/stir/SinogramIndices.inl new file mode 100644 index 000000000..bcb31b41f --- /dev/null +++ b/src/include/stir/SinogramIndices.inl @@ -0,0 +1,65 @@ +/*! + \file + \ingroup projdata + + \brief inline implementations for class stir::SinogramIndices + + \author Kris Thielemans + \author Sanida Mustafovic + \author PARAPET project + + +*/ +/* + Copyright (C) 2000 PARAPET partners + Copyright (C) 2000- 2009, Hammersmith Imanet Ltd + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license + + See STIR/LICENSE.txt for details +*/ + + +START_NAMESPACE_STIR + +SinogramIndices::SinogramIndices() +:SegmentIndices(),_axial_pos(0) + {} + +SinogramIndices::SinogramIndices( const int axial_pos_num,const int segment_num, const int timing_pos_num) + : SegmentIndices(segment_num, timing_pos_num),_axial_pos(axial_pos_num) + {} + +int +SinogramIndices::axial_pos_num() const +{ + return _axial_pos;} + + +int& +SinogramIndices::axial_pos_num() +{ return _axial_pos;} + +bool +SinogramIndices:: +operator<(const SinogramIndices& other) const +{ + return (_axial_pos< other._axial_pos) || + ((_axial_pos == other._axial_pos) && base_type::operator<(other)); +} + +bool +SinogramIndices:: +operator==(const SinogramIndices& other) const +{ + return (_axial_pos == other._axial_pos) && base_type::operator==(other); +} + +bool +SinogramIndices:: +operator!=(const SinogramIndices& other) const +{ + return !(*this == other); +} +END_NAMESPACE_STIR diff --git a/src/include/stir/Succeeded.h b/src/include/stir/Succeeded.h index 2b6c2d0da..57306aec6 100644 --- a/src/include/stir/Succeeded.h +++ b/src/include/stir/Succeeded.h @@ -16,6 +16,7 @@ /* Copyright (C) 2000 PARAPET partners Copyright (C) 2000- 2009, Hammersmith Imanet Ltd + Copyright (C) 2023, University College London This file is part of STIR. SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license @@ -35,15 +36,19 @@ START_NAMESPACE_STIR \code Succeeded f() { do_something; return Succeeded::yes; } void g() { if (f() == Succeeded::no) error("Error calling f"); } + // the latter can also be written as + void g2() { if (!f().succeeded()) error("Error calling f"); } \endcode */ class Succeeded { public: enum value { yes, no }; - Succeeded(const value& v) : v(v) {} + Succeeded(const value& v = yes) : v(v) {} bool operator==(const Succeeded &v2) const { return v == v2.v; } bool operator!=(const Succeeded &v2) const { return v != v2.v; } + //! convenience function returns if it is equal to Succeeded::yes + bool succeeded() const { return this->v == yes; } private: value v; }; diff --git a/src/include/stir/ViewSegmentNumbers.h b/src/include/stir/ViewSegmentNumbers.h index 67d89d942..2b359d0b2 100644 --- a/src/include/stir/ViewSegmentNumbers.h +++ b/src/include/stir/ViewSegmentNumbers.h @@ -4,78 +4,43 @@ \file \ingroup projdata - \brief Definition of class stir::ViewSegmentNumbers + \brief Definition of class stir::ViewSegmentNumbers, alias to stir::ViewgramIndices \author Kris Thielemans - \author Sanida Mustafovic - \author PARAPET project - + */ /* - Copyright (C) 2000 PARAPET partners - Copyright (C) 2000- 2009, Hammersmith Imanet Ltd + Copyright (C) 2023, University College London This file is part of STIR. - SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license + SPDX-License-Identifier: Apache-2.0 See STIR/LICENSE.txt for details */ - - #ifndef __stir_ViewSegmentNumbers_h__ #define __stir_ViewSegmentNumbers_h__ -#include "stir/common.h" +#include "stir/ViewgramIndices.h" START_NAMESPACE_STIR - +//! alias for ViewgramIndices /*! - \brief A very simple class to store view and segment numbers - \ingroup projdata + For backwards compatibility only. + + \deprecated */ -class ViewSegmentNumbers +class ViewSegmentNumbers : public ViewgramIndices { public: - - //! an empty constructor (sets everything to 0) - inline ViewSegmentNumbers(); - //! constructor taking view and segment number as arguments - inline ViewSegmentNumbers( const int view_num, const int segment_num, - const int tof_num = 0); - - //! get segment number for const objects - inline int segment_num() const; - //! get view number for const objects - inline int view_num() const; - //! get tof number for const objects - inline int tof_pos_num() const; - - //! get reference to segment number - inline int& segment_num(); - //! get reference to view number - inline int& view_num(); - //! get reference to timing position index - inline int& tof_pos_num(); - - - //! comparison operator, only useful for sorting - /*! order : (0,1) < (0,-1) < (1,1) ...*/ - inline bool operator<(const ViewSegmentNumbers& other) const; - - //! test for equality - inline bool operator==(const ViewSegmentNumbers& other) const; - inline bool operator!=(const ViewSegmentNumbers& other) const; - -private: - int segment; - int view; - int tof; - + using ViewgramIndices::ViewgramIndices; + // default constructor (needed for Visual Studio) + ViewSegmentNumbers() : ViewgramIndices() {} + ViewSegmentNumbers(const ViewgramIndices& ind) + : ViewgramIndices(ind) + {} }; END_NAMESPACE_STIR -#include "stir/ViewSegmentNumbers.inl" - #endif diff --git a/src/include/stir/Viewgram.h b/src/include/stir/Viewgram.h index 208e3b405..cbe870552 100644 --- a/src/include/stir/Viewgram.h +++ b/src/include/stir/Viewgram.h @@ -4,6 +4,7 @@ Copyright (C) 2000 PARAPET partners Copyright (C) 2000 - 2007-10-08, Hammersmith Imanet Ltd Copyright (C) 2011-07-01 - 2012, Kris Thielemans + Copyright (C) 2023, University College London This file is part of STIR. SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license @@ -28,7 +29,8 @@ #include "stir/Array.h" -#include "stir/ProjDataInfo.h" +#include "stir/ProjDataInfo.h" +#include "stir/ViewgramIndices.h" #include "stir/IndexRange.h" #include "stir/shared_ptr.h" @@ -39,7 +41,7 @@ START_NAMESPACE_STIR \ingroup projdata \brief A class for 2d projection data. - This represents a subset of the full projection. segment_num and view_num + This represents a subset of the full projection. SegmentIndices and view_num are fixed. */ @@ -60,15 +62,31 @@ class Viewgram : public Array<2,elemT> #endif public: + //! Construct from proj_data_info pointer and indices. Data are set to 0. + inline Viewgram(const shared_ptr& proj_data_info_ptr, + const ViewgramIndices& ind); + + //! Construct with data set to the array. + inline Viewgram(const Array<2,elemT>& p,const shared_ptr& proj_data_info_sptr, + const ViewgramIndices& ind); + //! Construct from proj_data_info pointer, view and segment number. Data are set to 0. + /*! + \deprecated Use version with ViewgramIndices instead + */ inline Viewgram(const shared_ptr& proj_data_info_ptr, const int v_num, const int s_num, const int t_num = 0); //! Construct with data set to the array. + /*! + \deprecated Use version with ViewgramIndices instead + */ inline Viewgram(const Array<2,elemT>& p,const shared_ptr& proj_data_info_ptr, const int v_num, const int s_num, const int t_num = 0); + //! Get indices + inline ViewgramIndices get_viewgram_indices() const; //! Get segment number inline int get_segment_num() const; //! Get number of views @@ -130,9 +148,7 @@ class Viewgram : public Array<2,elemT> private: shared_ptr proj_data_info_sptr; - int view_num; - int segment_num; - int timing_pos_num; + ViewgramIndices _indices; }; END_NAMESPACE_STIR diff --git a/src/include/stir/Viewgram.inl b/src/include/stir/Viewgram.inl index 02bbf4cf0..41d2d433a 100644 --- a/src/include/stir/Viewgram.inl +++ b/src/include/stir/Viewgram.inl @@ -16,6 +16,7 @@ /* Copyright (C) 2000 PARAPET partners Copyright (C) 2000- 2009, Hammersmith Imanet Ltd + Copyright (C) 2023, University College London This file is part of STIR. SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license @@ -27,20 +28,27 @@ START_NAMESPACE_STIR +template +ViewgramIndices +Viewgram::get_viewgram_indices() const +{ + return this->_indices; +} + template int Viewgram::get_segment_num() const -{ return segment_num; } +{ return this->_indices.segment_num(); } template int Viewgram::get_view_num() const -{ return view_num; } +{ return this->_indices.view_num(); } template int Viewgram::get_timing_pos_num() const -{ return timing_pos_num; } +{ return this->_indices.timing_pos_num(); } template int @@ -79,7 +87,7 @@ template Viewgram Viewgram::get_empty_copy(void) const { - Viewgram copy(proj_data_info_sptr, get_view_num(), get_segment_num(), get_timing_pos_num()); + Viewgram copy(proj_data_info_sptr, get_viewgram_indices()); return copy; } @@ -94,40 +102,54 @@ Viewgram::get_proj_data_info_sptr() const template Viewgram:: Viewgram(const Array<2,elemT>& p, - const shared_ptr& pdi_sptr, - const int v_num, const int s_num, const int t_num) + const shared_ptr& pdi_sptr, + const ViewgramIndices& ind) : Array<2,elemT>(p), proj_data_info_sptr(pdi_sptr), - view_num(v_num), segment_num(s_num), timing_pos_num(t_num) + _indices(ind) { - assert(view_num <= proj_data_info_sptr->get_max_view_num()); - assert(view_num >= proj_data_info_sptr->get_min_view_num()); + assert(ind.view_num() <= proj_data_info_sptr->get_max_view_num()); + assert(ind.view_num() >= proj_data_info_sptr->get_min_view_num()); // segment_num is already checked by doing get_max_axial_pos_num(s_num) - assert( get_min_axial_pos_num() == pdi_sptr->get_min_axial_pos_num(s_num)); - assert( get_max_axial_pos_num() == pdi_sptr->get_max_axial_pos_num(s_num)); + assert( get_min_axial_pos_num() == pdi_sptr->get_min_axial_pos_num(ind.segment_num())); + assert( get_max_axial_pos_num() == pdi_sptr->get_max_axial_pos_num(ind.segment_num())); assert( get_min_tangential_pos_num() == pdi_sptr->get_min_tangential_pos_num()); assert( get_max_tangential_pos_num() == pdi_sptr->get_max_tangential_pos_num()); } template Viewgram:: -Viewgram(const shared_ptr& pdi_sptr, - const int v_num, const int s_num, const int t_num) +Viewgram(const shared_ptr& pdi_sptr, + const ViewgramIndices& ind) : - Array<2,elemT>(IndexRange2D (pdi_sptr->get_min_axial_pos_num(s_num), - pdi_sptr->get_max_axial_pos_num(s_num), + Array<2,elemT>(IndexRange2D (pdi_sptr->get_min_axial_pos_num(ind.segment_num()), + pdi_sptr->get_max_axial_pos_num(ind.segment_num()), pdi_sptr->get_min_tangential_pos_num(), pdi_sptr->get_max_tangential_pos_num())), proj_data_info_sptr(pdi_sptr), - view_num(v_num), - segment_num(s_num), - timing_pos_num(t_num) + _indices(ind) { - assert(view_num <= proj_data_info_sptr->get_max_view_num()); - assert(view_num >= proj_data_info_sptr->get_min_view_num()); + assert(ind.view_num() <= proj_data_info_sptr->get_max_view_num()); + assert(ind.view_num() >= proj_data_info_sptr->get_min_view_num()); // segment_num is already checked by doing get_max_axial_pos_num(s_num) } +template +Viewgram:: +Viewgram(const Array<2,elemT>& p, + const shared_ptr& pdi_sptr, + const int v_num, const int s_num, const int t_num) + : + Viewgram(p, pdi_sptr, ViewgramIndices(v_num, s_num, t_num)) +{} + +template +Viewgram:: +Viewgram(const shared_ptr& pdi_sptr, + const int v_num, const int s_num, const int t_num) + : + Viewgram(pdi_sptr, ViewgramIndices(v_num, s_num, t_num)) +{} END_NAMESPACE_STIR diff --git a/src/include/stir/ViewgramIndices.h b/src/include/stir/ViewgramIndices.h new file mode 100644 index 000000000..6150ae22c --- /dev/null +++ b/src/include/stir/ViewgramIndices.h @@ -0,0 +1,65 @@ +// +// +/*! + \file + \ingroup projdata + + \brief Definition of class stir::ViewgramIndices + + \author Kris Thielemans + +*/ +/* + Copyright (C) 2023, University College London + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 + + See STIR/LICENSE.txt for details +*/ + +#ifndef __stir_ViewgramIndices_h__ +#define __stir_ViewgramIndices_h__ + +#include "stir/SegmentIndices.h" + +START_NAMESPACE_STIR + +/*! + \brief A very simple class to store all dincies to get a (2D) Viewgram + \ingroup projdata +*/ +class ViewgramIndices : public SegmentIndices +{ + typedef SegmentIndices base_type; + +public: + //! an empty constructor (sets everything to 0) + //TODOTOF remove default + inline ViewgramIndices(); + //! constructor taking view and segment number as arguments + inline ViewgramIndices(const int view_num, const int segment_num, const int timing_pos_num=0); + + //! get view number for const objects + inline int view_num() const; + + //! get reference to view number + inline int& view_num(); + + //! comparison operator, only useful for sorting + /*! order : (0,1) < (0,-1) < (1,1) ...*/ + inline bool operator<(const ViewgramIndices& other) const; + + //! test for equality + inline bool operator==(const ViewgramIndices& other) const; + inline bool operator!=(const ViewgramIndices& other) const; + +private: + int _view; +}; + +END_NAMESPACE_STIR + +#include "stir/ViewgramIndices.inl" + +#endif diff --git a/src/include/stir/ViewgramIndices.inl b/src/include/stir/ViewgramIndices.inl new file mode 100644 index 000000000..b6bc88c65 --- /dev/null +++ b/src/include/stir/ViewgramIndices.inl @@ -0,0 +1,64 @@ +/*! + \file + \ingroup projdata + + \brief inline implementations for class stir::ViewgramIndices + + \author Kris Thielemans + \author Sanida Mustafovic + \author PARAPET project + + +*/ +/* + Copyright (C) 2000 PARAPET partners + Copyright (C) 2000- 2009, Hammersmith Imanet Ltd + This file is part of STIR. + + SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license + + See STIR/LICENSE.txt for details +*/ + +START_NAMESPACE_STIR + +ViewgramIndices::ViewgramIndices() + : SegmentIndices(), + _view(0) +{} + +ViewgramIndices::ViewgramIndices(const int view_num, const int segment_num, const int timing_pos_num) + : SegmentIndices(segment_num, timing_pos_num), + _view(view_num) +{} + +int +ViewgramIndices::view_num() const +{ + return _view; +} + +int& +ViewgramIndices::view_num() +{ + return _view; +} + +bool +ViewgramIndices::operator<(const ViewgramIndices& other) const +{ + return (_view < other._view) || ((_view == other._view) && base_type::operator<(other)); +} + +bool +ViewgramIndices::operator==(const ViewgramIndices& other) const +{ + return (_view == other._view) && base_type::operator==(other); +} + +bool +ViewgramIndices::operator!=(const ViewgramIndices& other) const +{ + return !(*this == other); +} +END_NAMESPACE_STIR diff --git a/src/include/stir/common.h b/src/include/stir/common.h index 3f2315aee..9ada91aa8 100644 --- a/src/include/stir/common.h +++ b/src/include/stir/common.h @@ -34,34 +34,23 @@
      • macros for namespace support: - #defines STIR_NO_NAMESPACES if the compiler does not support namespaces - #defines START_NAMESPACE_STIR etc. - -
      • #defines STIR_NO_COVARIANT_RETURN_TYPES when the compiler does not - support virtual functions of a derived class differing only in the return - type. - Define when your compiler does not handle the following: - \code - class A { virtual A* f();} - class B:A { virtual B* f(); } - \endcode + \c \#defines \c START_NAMESPACE_STIR etc. +
      • preprocessor definitions which attempt to determine the operating system this is going to run on. - use as #ifdef __OS_WIN__ ... #elif ... #endif + use as \c "#ifdef __OS_WIN__ ... #elif ... #endif" Possible values are __OS_WIN__, __OS_MAC__, __OS_VAX__, __OS_UNIX__ The __OS_UNIX__ case has 'subbranches'. At the moment we attempt to find out on __OS_AIX__, __OS_SUN__, __OS_OSF__, __OS_LINUX__. (If the attempts fail to determine the correct OS, you can pass the correct value as a preprocessor definition to the compiler) - -
      • #includes cstdio, cstdlib, cstring, cmath - -
      • a trick to get ANSI C++ 'for' scoping rules work for compilers - which do not follow the new standard - -
      • #ifdef STIR_ASSERT, then define our own assert, else include <cassert> +
        • +
      • \c \#includes cstdio, cstdlib, cstring, cmath +
        • +
      • \c \#ifdef \c STIR_ASSERT, then define our own assert, else include <cassert> +

      Speeding up std::copy

      @@ -69,6 +58,7 @@
      • For old compilers (check the source!), overloads of std::copy for built-in types to use memmove (so it's faster) +

      stir namespace members declared here

      @@ -83,8 +73,6 @@

      stir include files included here

      • stir/config.h sets various preprocessor defines (generated from STIRConfig.in)
        • -
      • stir/error.h defining stir::error
        • -
      • stir/warning.h defining stir::warning
      */ #include "stir/config.h" @@ -96,22 +84,12 @@ #include //*************** namespace macros -#ifndef STIR_NO_NAMESPACES - # define START_NAMESPACE_STIR namespace stir { # define END_NAMESPACE_STIR } # define USING_NAMESPACE_STIR using namespace stir; # define START_NAMESPACE_STD namespace std { # define END_NAMESPACE_STD } # define USING_NAMESPACE_STD using namespace std; -#else -# define START_NAMESPACE_STIR -# define END_NAMESPACE_STIR -# define USING_NAMESPACE_STIR -# define START_NAMESPACE_STD -# define END_NAMESPACE_STD -# define USING_NAMESPACE_STD -#endif //*************** define __OS_xxx__ @@ -239,23 +217,6 @@ END_NAMESPACE_STD #endif // #ifdef STIR_SPEED_UP_STD_COPY -//*************** for() scope trick - -/* ANSI C++ (re)defines the scope of variables declared in a for() statement. - Example: the 'i' variable has scope only within the for statement. - for (int i=0; ...) - do_something; - The next trick (by AZ) solves this problem. - At the moment, we only need it for VC++ - */ - -#if defined(STIR_ENABLE_FOR_SCOPE_WORKAROUND) -# ifndef for -# define for if (0) ; else for -# else -# error 'for' is already #defined to something -# endif -#endif //*************** assert diff --git a/src/include/stir/recon_buildblock/BinNormalisationFromECAT8.h b/src/include/stir/recon_buildblock/BinNormalisationFromECAT8.h index 1bc57505c..d5c894a64 100644 --- a/src/include/stir/recon_buildblock/BinNormalisationFromECAT8.h +++ b/src/include/stir/recon_buildblock/BinNormalisationFromECAT8.h @@ -1,6 +1,6 @@ /* Copyright (C) 2000-2007, Hammersmith Imanet Ltd - Copyright (C) 2013-2014 University College London + Copyright (C) 2013-2014, 2020, 2023 University College London Largely a copy of the ECAT7 version. @@ -49,20 +49,40 @@ START_NAMESPACE_ECAT \par Parsing example \verbatim - Bin Normalisation type := from ecat7 + Bin Normalisation type := from ecat8 Bin Normalisation From ECAT8:= - normalisation filename:= myfile.hn + normalisation filename:= myfile.n.hdr ; next keywords can be used to switch off some of the normalisation components - ; do not use unless you know why + ; do not use unless you know why. + ; Default values are indicated below (i.e. use all of them) ; use_gaps:=1 ; use_detector_efficiencies:=1 ; use_dead_time:=1 ; use_geometric_factors:=1 ; use_crystal_interference_factors:=1 + ; use_axial_effects_factors:=1 + + ; keyword that can be used to write the components to a separate text files for debugging + ; files are written in the current directory and are called geom_out.txt etc. + ; write_components_to_file := 0 End Bin Normalisation From ECAT8:= \endverbatim + \par More information + + Siemens stores `axial effects`, i.e. one number per sinogram. This normally limits the use of the + file to data that have been acquired with the same `span`, which is usually 11 for present Siemens scanners. + + We work around this in 2 ways: + - we assume that the `axial_effects_factor` is the same for every ring pair contributing to a particular + Siemens sinogram + - if there is no corresponding Siemens sinogram (i.e. the norm file has been acquired with a particular + maximum ring difference, smaller than what is actually possible with the scanner), we use an + `axial_effects_factor` of 1. This should be reasonable as the numbers are around 1 (on the mMR). + + This strategy allows us to give normalisation factor for span=1 data, even if the norm file is for span=11. + \todo dead-time is not yet implemented @@ -94,8 +114,9 @@ class BinNormalisationFromECAT8 : bool use_dead_time() const; bool use_geometric_factors() const; bool use_crystal_interference_factors() const; + bool use_axial_effects_factors() const; -private: + private: Array<1,float> axial_t1_array; Array<1,float> axial_t2_array; Array<1,float> trans_t1_array; @@ -103,14 +124,20 @@ class BinNormalisationFromECAT8 : Array<2,float> geometric_factors; Array<2,float> efficiency_factors; Array<2,float> crystal_interference_factors; + Array<1,float> axial_effects; + //! lookup table from STIR ring-pair to a Siemens sinogram-index + Array<2,int> sino_index; + //! number of sinograms in Siemens sinogram (span=11?) + int num_Siemens_sinograms; + shared_ptr scanner_ptr; int num_transaxial_crystals_per_block; // TODO move to Scanner int num_axial_blocks_per_singles_unit; shared_ptr proj_data_info_ptr; + shared_ptr norm_proj_data_info_sptr; ProjDataInfoCylindricalNoArcCorr const * proj_data_info_cyl_ptr; shared_ptr proj_data_info_cyl_uncompressed_ptr; - int span; int mash; int num_blocks_per_singles_unit; float calib_factor, cross_calib_factor; @@ -120,11 +147,16 @@ class BinNormalisationFromECAT8 : bool _use_dead_time; bool _use_geometric_factors; bool _use_crystal_interference_factors; + bool _use_axial_effects_factors; + bool _write_components_to_file; void read_norm_data(const string& filename); float get_dead_time_efficiency ( const DetectionPosition<>& det_pos, const double start_time, const double end_time) const; + //! initialise sino_index and num_Siemens_sinograms + void construct_sino_lookup_table(); + float find_axial_effects(int ring1, int ring2) const; // parsing stuff virtual void set_defaults() override; virtual void initialise_keymap() override; diff --git a/src/include/stir/recon_buildblock/ForwardProjectorByBinUsingRayTracing.h b/src/include/stir/recon_buildblock/ForwardProjectorByBinUsingRayTracing.h index 057ce8abc..bc07b2e94 100644 --- a/src/include/stir/recon_buildblock/ForwardProjectorByBinUsingRayTracing.h +++ b/src/include/stir/recon_buildblock/ForwardProjectorByBinUsingRayTracing.h @@ -44,8 +44,8 @@ class ProjDataInfoCylindrical; \brief This class implements forward projection using Siddon's algorithm for ray tracing. That is, it computes length of intersection with the voxels. - Currently, the LOIs are divided by voxel_size.x(), unless NEWSCALE is - #defined during compilation time of ForwardProjectorByBinUsingRayTracing_Siddon.cxx. + Currently, the LOIs are divided by voxel_size.x(), unless \c NEWSCALE is + \c \#defined during compilation time of ForwardProjectorByBinUsingRayTracing_Siddon.cxx. If the z voxel size is exactly twice the sampling in axial direction, multiple LORs are used, to avoid missing voxels. (TODOdoc describe how). diff --git a/src/include/stir/recon_buildblock/GeneralisedObjectiveFunction.h b/src/include/stir/recon_buildblock/GeneralisedObjectiveFunction.h index b6fbcad91..021089257 100644 --- a/src/include/stir/recon_buildblock/GeneralisedObjectiveFunction.h +++ b/src/include/stir/recon_buildblock/GeneralisedObjectiveFunction.h @@ -265,9 +265,6 @@ class GeneralisedObjectiveFunction: return exam_info_uptr; } - //! Return the status of TOF - bool get_tof_status() const; - //! Attempts to change the number of subsets. /*! \return The number of subsets that will be used later, which is not guaranteed to be what you asked for. */ @@ -335,9 +332,6 @@ class GeneralisedObjectiveFunction: protected: int num_subsets; - //! If set TOF information will be taken into account. - bool use_tof; - shared_ptr > prior_sptr; //! sets any default values diff --git a/src/include/stir/recon_buildblock/PLSPrior.h b/src/include/stir/recon_buildblock/PLSPrior.h index f2e363991..018d47c44 100644 --- a/src/include/stir/recon_buildblock/PLSPrior.h +++ b/src/include/stir/recon_buildblock/PLSPrior.h @@ -48,7 +48,7 @@ START_NAMESPACE_STIR \phi(f) = \sqrt{\alpha^2 + |\nabla f|^2 - {\langle\nabla f,\xi\rangle}^2} \f] where \f$ f \f$ is the PET image, - \f$\nabla \f$ is the finite difference operator (\textrm{not} taking voxel-sizes into account) and + \f$\nabla \f$ is the finite difference operator (\b not taking voxel-sizes into account) and \f$ \xi \f$ is the normalised gradient of the anatomical image calculated as follows: \f[ diff --git a/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.h b/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.h index 74b094b18..2b7733d66 100644 --- a/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.h +++ b/src/include/stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.h @@ -334,7 +334,7 @@ public RegisteredParsingObject additive_proj_data_sptr; shared_ptr normalisation_sptr; - + // TODO doc int frame_num; std::string frame_definition_filename; @@ -390,6 +390,18 @@ public RegisteredParsingObject sens_proj_data_info_sptr; + //! flag to check if we set-up normalisation with the original data or not + mutable bool latest_setup_norm_was_with_orig_data; + //! flag to check if normalisation has been set_up already or not + mutable bool norm_already_setup = false; + + //! helper function to makre sure we set-up normalisation_sptr correctly + void ensure_norm_is_set_up(bool for_original_data = true) const; + //! convenience for ensure_norm_is_set_up(false) + void ensure_norm_is_set_up_for_sensitivity() const; //!@} #if 0 void diff --git a/src/include/stir/recon_buildblock/ProjMatrixByBinPinholeSPECTUB.h b/src/include/stir/recon_buildblock/ProjMatrixByBinPinholeSPECTUB.h index 13305f625..cf9dacd99 100644 --- a/src/include/stir/recon_buildblock/ProjMatrixByBinPinholeSPECTUB.h +++ b/src/include/stir/recon_buildblock/ProjMatrixByBinPinholeSPECTUB.h @@ -3,17 +3,7 @@ Copyright (C) 2021, University College London This file is part of STIR. - Licensed under the Apache License, Version 2.0 (the "License"); - you may not use this file except in compliance with the License. - You may obtain a copy of the License at - - http://www.apache.org/licenses/LICENSE-2.0 - - Unless required by applicable law or agreed to in writing, software - distributed under the License is distributed on an "AS IS" BASIS, - WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - See the License for the specific language governing permissions and - limitations under the License. + SPDX-License-Identifier: Apache-2.0 See STIR/LICENSE.txt for details */ @@ -53,7 +43,7 @@ class Bin; \warning this class currently only works with VoxelsOnCartesianGrid. - \Sample parameter file + \par Sample parameter file \verbatim Projection Matrix By Bin Pinhole SPECT UB Parameters:= @@ -83,7 +73,7 @@ class Bin; End Projection Matrix By Bin Pinhole SPECT UB Parameters:= \endverbatim - \Sample detector file + \par Sample detector file \verbatim Information of detector @@ -104,7 +94,7 @@ class Bin; \#…………until here……………\# \endverbatim - \Sample collimator file + \par Sample collimator file \verbatim Information of collimator Comments are allowed here or anywhere in lines not containing parameters. diff --git a/src/include/stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h b/src/include/stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h index 14609deb4..a81577a72 100644 --- a/src/include/stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h +++ b/src/include/stir/recon_buildblock/ProjMatrixByBinUsingRayTracing.h @@ -38,8 +38,8 @@ class ProjDataInfo; \brief Computes projection matrix elements for VoxelsOnCartesianGrid images by using a Length of Intersection (LOI) model. - Currently, the LOIs are divided by voxel_size.x(), unless NEWSCALE is - #defined during compilation time of ProjMatrixByBinUsingRayTracing.cxx. + Currently, the LOIs are divided by voxel_size.x(), unless \c NEWSCALE is + \c \#defined during compilation time of ProjMatrixByBinUsingRayTracing.cxx. It is possible to use multiple LORs in tangential direction. The result will then be the average of the various contributions. Currently all these diff --git a/src/include/stir_experimental/modelling/BloodFrameData.inl b/src/include/stir_experimental/modelling/BloodFrameData.inl index e526910fe..5fd4a1a1d 100644 --- a/src/include/stir_experimental/modelling/BloodFrameData.inl +++ b/src/include/stir_experimental/modelling/BloodFrameData.inl @@ -25,7 +25,8 @@ START_NAMESPACE_STIR BloodFrameData::BloodFrameData() { } -//! constructor giving a vector //ChT::ToDO: Better to use iterators +//! constructor giving a vector +//ChT::ToDO: Better to use iterators BloodFrameData::BloodFrameData(const std::vector & blood_plot) {this->_blood_plot=blood_plot;} diff --git a/src/include/stir_experimental/motion/RigidObject3DTransformation.h b/src/include/stir_experimental/motion/RigidObject3DTransformation.h index c3432a83b..7fc0aa723 100644 --- a/src/include/stir_experimental/motion/RigidObject3DTransformation.h +++ b/src/include/stir_experimental/motion/RigidObject3DTransformation.h @@ -143,7 +143,7 @@ class RigidObject3DTransformation //! Transform bin from some projection data /*! Finds 'closest' (in some sense) bin to the transformed LOR. - if NEW_ROT is not #defined at compilation time, + if \c NEW_ROT is not \c \#defined at compilation time, it will throw an exception when arc-corrected data is used.*/ void transform_bin(Bin& bin,const ProjDataInfo& out_proj_data_info, const ProjDataInfo& in_proj_data_info) const; diff --git a/src/include/stir_experimental/recon_buildblock/ProjMatrixByDenselOnCartesianGridUsingElement.h b/src/include/stir_experimental/recon_buildblock/ProjMatrixByDenselOnCartesianGridUsingElement.h index 47baf865d..a1b8304b6 100644 --- a/src/include/stir_experimental/recon_buildblock/ProjMatrixByDenselOnCartesianGridUsingElement.h +++ b/src/include/stir_experimental/recon_buildblock/ProjMatrixByDenselOnCartesianGridUsingElement.h @@ -32,8 +32,8 @@ template class DiscretisedDensity; \brief Computes projection matrix elements for VoxelsOnCartesianGrid images by using a Length of Intersection (LOI) model. - Currently, the LOIs are divided by voxel_size.x(), unless NEWSCALE is - #defined during compilation time of ProjMatrixByDenselOnCartesianGridUsingElement.cxx. + Currently, the LOIs are divided by voxel_size.x(), unless \c NEWSCALE is + \c \#defined during compilation time of ProjMatrixByDenselOnCartesianGridUsingElement.cxx. If the z voxel size is exactly twice the sampling in axial direction, multiple LORs are used, to avoid missing voxels. (TODOdoc describe how). diff --git a/src/include/stir_experimental/recon_buildblock/ProjMatrixByDenselUsingRayTracing.h b/src/include/stir_experimental/recon_buildblock/ProjMatrixByDenselUsingRayTracing.h index e3bb0adef..a1da9e19e 100644 --- a/src/include/stir_experimental/recon_buildblock/ProjMatrixByDenselUsingRayTracing.h +++ b/src/include/stir_experimental/recon_buildblock/ProjMatrixByDenselUsingRayTracing.h @@ -32,8 +32,8 @@ class DataSymmetriesForDensels_PET_CartesianGrid; \brief Computes projection matrix elements for VoxelsOnCartesianGrid images by using a Length of Intersection (LOI) model. - Currently, the LOIs are divided by voxel_size.x(), unless NEWSCALE is - #defined during compilation time of ProjMatrixByDenselUsingRayTracing.cxx. + Currently, the LOIs are divided by voxel_size.x(), unless \c NEWSCALE is + \c \#defined during compilation time of ProjMatrixByDenselUsingRayTracing.cxx. If the z voxel size is exactly twice the sampling in axial direction, multiple LORs are used, to avoid missing voxels. (TODOdoc describe how). diff --git a/src/iterative/KOSMAPOSL/CMakeLists.txt b/src/iterative/KOSMAPOSL/CMakeLists.txt index 3858baee1..cc59bb67f 100644 --- a/src/iterative/KOSMAPOSL/CMakeLists.txt +++ b/src/iterative/KOSMAPOSL/CMakeLists.txt @@ -23,5 +23,5 @@ if (STIR_MPI) SET_PROPERTY(TARGET KOSMAPOSL PROPERTY LINK_FLAGS ${MPI_CXX_LINK_FLAGS}) endif() -target_link_libraries(iterative_KOSMAPOSL recon_buildblock ) +target_link_libraries(iterative_KOSMAPOSL PUBLIC recon_buildblock ) diff --git a/src/iterative/OSMAPOSL/CMakeLists.txt b/src/iterative/OSMAPOSL/CMakeLists.txt index 64ef89a2e..d1f26ec13 100644 --- a/src/iterative/OSMAPOSL/CMakeLists.txt +++ b/src/iterative/OSMAPOSL/CMakeLists.txt @@ -23,5 +23,5 @@ if (STIR_MPI) SET_PROPERTY(TARGET OSMAPOSL PROPERTY LINK_FLAGS ${MPI_CXX_LINK_FLAGS}) endif() -target_link_libraries(iterative_OSMAPOSL recon_buildblock ) +target_link_libraries(iterative_OSMAPOSL PUBLIC recon_buildblock ) diff --git a/src/iterative/OSSPS/CMakeLists.txt b/src/iterative/OSSPS/CMakeLists.txt index d5f7aafdf..9e7f77a05 100644 --- a/src/iterative/OSSPS/CMakeLists.txt +++ b/src/iterative/OSSPS/CMakeLists.txt @@ -24,4 +24,4 @@ if (STIR_MPI) SET_PROPERTY(TARGET OSSPS PROPERTY LINK_FLAGS ${MPI_CXX_LINK_FLAGS}) endif() -target_link_libraries(iterative_OSSPS recon_buildblock ) +target_link_libraries(iterative_OSSPS PUBLIC recon_buildblock ) diff --git a/src/listmode_buildblock/CMakeLists.txt b/src/listmode_buildblock/CMakeLists.txt index b21d21d3d..baed95da4 100644 --- a/src/listmode_buildblock/CMakeLists.txt +++ b/src/listmode_buildblock/CMakeLists.txt @@ -57,4 +57,4 @@ endif() include(stir_lib_target) -target_link_libraries(listmode_buildblock data_buildblock ) +target_link_libraries(listmode_buildblock PUBLIC data_buildblock ) diff --git a/src/modelling_buildblock/CMakeLists.txt b/src/modelling_buildblock/CMakeLists.txt index aab6d86d8..aed39c0c1 100644 --- a/src/modelling_buildblock/CMakeLists.txt +++ b/src/modelling_buildblock/CMakeLists.txt @@ -15,4 +15,4 @@ set(${dir_LIB_SOURCES} include(stir_lib_target) -target_link_libraries(modelling_buildblock buildblock IO) +target_link_libraries(modelling_buildblock PUBLIC buildblock IO) diff --git a/src/numerics_buildblock/CMakeLists.txt b/src/numerics_buildblock/CMakeLists.txt index 823e78daf..c6e8f4158 100644 --- a/src/numerics_buildblock/CMakeLists.txt +++ b/src/numerics_buildblock/CMakeLists.txt @@ -13,4 +13,4 @@ set(${dir_LIB_SOURCES} include(stir_lib_target) -target_link_libraries(${dir} buildblock) +target_link_libraries(${dir} PUBLIC buildblock) diff --git a/src/recon_buildblock/BackProjectorByBin.cxx b/src/recon_buildblock/BackProjectorByBin.cxx index 1ab6750cb..59a35bb22 100644 --- a/src/recon_buildblock/BackProjectorByBin.cxx +++ b/src/recon_buildblock/BackProjectorByBin.cxx @@ -166,6 +166,8 @@ back_project(DiscretisedDensity<3,float>& density, { ViewSegmentNumbers vs(iter->get_view_num(), iter->get_segment_num()); get_symmetries_used()->find_basic_view_segment_numbers(vs); + // TODOTOF find_basic_view_segment_numbers doesn't fill in timing_pos_num + vs.timing_pos_num() = basic_vs.timing_pos_num(); if (vs != basic_vs) error("BackProjectorByBin::back_project called with incorrect related_viewgrams. Problem with symmetries!\n"); } @@ -217,12 +219,13 @@ BackProjectorByBin::back_project(const ProjData& proj_data, int subset_num, int RelatedViewgrams viewgrams; #pragma omp critical (BACKPROJECTORBYBIN_GETVIEWGRAMS) viewgrams = proj_data.get_related_viewgrams(vs, symmetries_sptr, false, k); + info(boost::format("Processing view %1% of segment %2%, TOF bin %3%") % vs.view_num() % vs.segment_num() % k, 3); #else const RelatedViewgrams viewgrams = proj_data.get_related_viewgrams(vs, symmetries_sptr, false, k); + info(boost::format("Processing view %1% of segment %2%, TOF bin %3%") % vs.view_num() % vs.segment_num() % k, 3); #endif - info(boost::format("Processing view %1% of segment %2%") % vs.view_num() % vs.segment_num(), 2); back_project(viewgrams); } } @@ -285,6 +288,8 @@ back_project(const RelatedViewgrams& viewgrams, { ViewSegmentNumbers vs(iter->get_view_num(), iter->get_segment_num()); get_symmetries_used()->find_basic_view_segment_numbers(vs); + // TODOTOF find_basic_view_segment_numbers doesn't fill in timing_pos_num + vs.timing_pos_num() = basic_vs.timing_pos_num(); if (vs != basic_vs) error("BackProjectorByBin::back_project called with incorrect related_viewgrams. Problem with symmetries!\n"); } diff --git a/src/recon_buildblock/BackProjectorByBinUsingInterpolation_linear.cxx b/src/recon_buildblock/BackProjectorByBinUsingInterpolation_linear.cxx index 14782cb02..9061e5979 100644 --- a/src/recon_buildblock/BackProjectorByBinUsingInterpolation_linear.cxx +++ b/src/recon_buildblock/BackProjectorByBinUsingInterpolation_linear.cxx @@ -9,7 +9,7 @@ stir::BackProjectorByBinUsingInterpolation, for the case of piecewise linear interpolation. - \warning This #includes BackProjectorByBinUsingInterpolation_3DCho.cxx + \warning This \c \#includes BackProjectorByBinUsingInterpolation_3DCho.cxx This very ugly way of including a .cxx file is used to avoid replication of a lot of (difficult) code. diff --git a/src/recon_buildblock/BackProjectorByBinUsingInterpolation_piecewise_linear.cxx b/src/recon_buildblock/BackProjectorByBinUsingInterpolation_piecewise_linear.cxx index 716d1cef8..aa5ed32f4 100644 --- a/src/recon_buildblock/BackProjectorByBinUsingInterpolation_piecewise_linear.cxx +++ b/src/recon_buildblock/BackProjectorByBinUsingInterpolation_piecewise_linear.cxx @@ -9,7 +9,7 @@ stir::BackProjectorByBinUsingInterpolation, for the case of piecewise linear interpolation. - \warning This #includes BackProjectorByBinUsingInterpolation_3DCho.cxx + \warning This \c \#includes BackProjectorByBinUsingInterpolation_3DCho.cxx This very ugly way of including a .cxx file is used to avoid replication of a lot of (difficult) code. diff --git a/src/recon_buildblock/BinNormalisation.cxx b/src/recon_buildblock/BinNormalisation.cxx index 776401558..3f6919b9c 100644 --- a/src/recon_buildblock/BinNormalisation.cxx +++ b/src/recon_buildblock/BinNormalisation.cxx @@ -61,11 +61,11 @@ get_exam_info_sptr() const Succeeded BinNormalisation:: -set_up(const shared_ptr& exam_info_sptr, const shared_ptr& proj_data_info_sptr_v ) +set_up(const shared_ptr& exam_info_sptr_v, const shared_ptr& proj_data_info_sptr_v ) { _already_set_up = true; - this->proj_data_info_sptr = proj_data_info_sptr_v->create_shared_clone(); - this->exam_info_sptr=exam_info_sptr; + this->proj_data_info_sptr = proj_data_info_sptr_v; + this->exam_info_sptr = exam_info_sptr_v; return Succeeded::yes; } diff --git a/src/recon_buildblock/BinNormalisationFromECAT8.cxx b/src/recon_buildblock/BinNormalisationFromECAT8.cxx index 613d170cc..4c6ed5811 100644 --- a/src/recon_buildblock/BinNormalisationFromECAT8.cxx +++ b/src/recon_buildblock/BinNormalisationFromECAT8.cxx @@ -1,6 +1,7 @@ /* Copyright (C) 2002-2011, Hammersmith Imanet Ltd - Copyright (C) 2013-2014, 2019 University College London + Copyright (C) 2013-2014, 2019, 2020, 2021 University College London + Copyright (C) 2020, National Physical Laboratory This file contains is based on information supplied by Siemens but is distributed with their consent. @@ -22,6 +23,7 @@ \author Kris Thielemans \author Sanida Mustafovic + \author Daniel Deidda */ @@ -36,9 +38,10 @@ #include "stir/Bin.h" #include "stir/display.h" #include "stir/IO/read_data.h" -#include "stir/IO/InterfileHeader.h" +#include "stir/IO/InterfileHeaderSiemens.h" #include "stir/ByteOrder.h" #include "stir/is_null_ptr.h" +#include "stir/utilities.h" #include "stir/warning.h" #include "stir/error.h" #include @@ -160,7 +163,9 @@ BinNormalisationFromECAT8::set_defaults() this->_use_detector_efficiencies = true; this->_use_dead_time = false; this->_use_geometric_factors = true; - this->_use_crystal_interference_factors = true; + this->_use_crystal_interference_factors = true; + this->_use_axial_effects_factors = true; + this->_write_components_to_file = false; } void @@ -178,6 +183,8 @@ initialise_keymap() //this->parser.add_key("use_dead_time", &this->_use_dead_time); this->parser.add_key("use_geometric_factors", &this->_use_geometric_factors); this->parser.add_key("use_crystal_interference_factors", &this->_use_crystal_interference_factors); + this->parser.add_key("use_axial_effects_factors", &this->_use_axial_effects_factors); + this->parser.add_key("write_components_to_file", &this->_write_components_to_file); this->parser.add_stop_key("End Bin Normalisation From ECAT8"); } @@ -229,12 +236,20 @@ set_up(const shared_ptr &exam_info_sptr_v, const shared_ptrget_max_ring_difference(0) - - proj_data_info_cyl_ptr->get_min_ring_difference(0) + 1; - // TODO insert check all other segments are the same + if (this->use_axial_effects_factors()) + { + const int data_max_ring_diff = + proj_data_info_cyl_ptr->get_max_ring_difference(proj_data_info_cyl_ptr->get_max_segment_num()); + auto norm_proj_data_info_no_arccorr_ptr = + dynamic_cast(norm_proj_data_info_sptr.get()); + const int norm_max_ring_diff = + norm_proj_data_info_no_arccorr_ptr->get_max_ring_difference(norm_proj_data_info_no_arccorr_ptr->get_max_segment_num()); + if (data_max_ring_diff > norm_max_ring_diff) + warning("ECAT8 norm axial effects given only for ring diff up to " + std::to_string(norm_max_ring_diff) + + ". I will use 1 for larger ring difference."); + } - mash = scanner_ptr->get_num_detectors_per_ring()/2/proj_data_info_ptr->get_num_views(); + this->mash = scanner_ptr->get_num_detectors_per_ring()/2/proj_data_info_ptr->get_num_views(); return Succeeded::yes; } @@ -244,65 +259,17 @@ BinNormalisationFromECAT8:: read_norm_data(const string& filename) { -#if 0 -MatrixFile* mptr = matrix_open(filename.c_str(), MAT_READ_ONLY, Norm3d); - if (mptr == 0) - error("BinNormalisationFromECAT8: error opening %s\n", filename.c_str()); + InterfileNormHeaderSiemens norm_parser; + norm_parser.parse(filename.c_str()); - scanner_ptr.reset( - find_scanner_from_ECAT_system_type(mptr->mhptr->system_type)); - - MatrixData* matrix = matrix_read( mptr, mat_numcod (1, 1, 1, 0, 0), - Norm3d /*= read data as well */); - - num_transaxial_crystals_per_block = nrm_subheader_ptr->num_transaxial_crystals ; -#endif -#if 0 - InterfileHeader interfile_parser; - ignore_key("data format"); - interfile_parser.parse(filename.c_str()); - -#else - KeyParser parser; - std::string originating_system; - std::string data_file_name; - int num_buckets; - { - parser.add_start_key("INTERFILE"); - parser.add_stop_key("END OF INTERFILE"); // add this for safety (even though it isn't always there) - parser.add_key("originating_system", &originating_system); - parser.add_key("name_of_data_file", &data_file_name); - parser.add_key("%number of buckets", &num_buckets); - parser.add_key("%scanner quantification factor (Bq*s/ECAT counts)",& calib_factor); - parser.add_key("%cross calibration factor",& cross_calib_factor); - parser.parse(filename.c_str()); - } -#endif - // remove trailing \r - std::string s=/*interfile_parser.*/originating_system; - s.erase( std::remove_if( s.begin(), s.end(), isspace ), s.end() ); - /*interfile_parser.*/originating_system=s; - s=/*interfile_parser.*/data_file_name; - s.erase( std::remove_if( s.begin(), s.end(), isspace ), s.end() ); - /*interfile_parser.*/data_file_name=s; - - this->scanner_ptr.reset(Scanner::get_scanner_from_name(/*interfile_parser.*/originating_system)); - switch(this->scanner_ptr->get_type()) - { - //case Scanner::E1080: - case Scanner::Siemens_mCT: - case Scanner::Siemens_mMR: - break; - default: - error(boost::format("Unknown originating_system '%s', when parsing file '%s'") % /*interfile_parser.*/originating_system % filename ); - } + this->norm_proj_data_info_sptr = norm_parser.data_info_ptr; + this->scanner_ptr = norm_parser.data_info_ptr->get_scanner_sptr(); - char directory_name[max_filename_length]; + char directory_name[max_filename_length]; get_directory_name(directory_name, filename.c_str()); char full_data_file_name[max_filename_length]; - strcpy(full_data_file_name, data_file_name.c_str()); + strcpy(full_data_file_name, norm_parser.data_file_name.c_str()); prepend_directory_name(full_data_file_name, directory_name); - num_transaxial_crystals_per_block = scanner_ptr->get_num_transaxial_crystals_per_block(); // Calculate the number of axial blocks per singles unit and // total number of blocks per singles unit. @@ -348,7 +315,8 @@ MatrixFile* mptr = matrix_open(filename.c_str(), MAT_READ_ONLY, Norm3d); /*num_tangential_poss=*/scanner_ptr->get_max_num_non_arccorrected_bins(), //XXXnrm_subheader_ptr->num_r_elements, /*arc_corrected =*/false) )); - + + this->construct_sino_lookup_table(); /* Extract geometrical & crystal interference, and crystal efficiencies from the normalisation data. @@ -377,21 +345,26 @@ MatrixFile* mptr = matrix_open(filename.c_str(), MAT_READ_ONLY, Norm3d); efficiency_factors = Array<2,float>(IndexRange2D(0,scanner_ptr->get_num_rings()-1, 0, scanner_ptr->get_num_detectors_per_ring()-1)); - -#if 0 - int geom_test = nrm_subheader_ptr->num_geo_corr_planes * (max_tang_pos_num-min_tang_pos_num +1); - int cry_inter = num_transaxial_crystals_per_block * (max_tang_pos_num-min_tang_pos_num +1); - int eff_test = scanner_ptr->get_num_detectors_per_ring() * scanner_ptr->get_num_rings(); -#endif + axial_effects = + Array<1,float>(num_Siemens_sinograms); - std::ifstream binary_data(full_data_file_name, std::ios::binary | std::ios::in); + std::ifstream binary_data; + open_read_binary(binary_data, full_data_file_name); if (read_data(binary_data, geometric_factors, ByteOrder::little_endian) != Succeeded::yes) error("failed reading geo factors from '%s'", full_data_file_name); + if (binary_data.tellg() != std::streampos(norm_parser.data_offset_each_dataset[1])) + error("Error reading ECAT8 norm file: wrong offset after component 1"); if (read_data(binary_data, crystal_interference_factors, ByteOrder::little_endian) != Succeeded::yes) error("failed reading crystal_interference_factors from '%s'", full_data_file_name); + if (binary_data.tellg() != std::streampos(norm_parser.data_offset_each_dataset[2])) + error("Error reading ECAT8 norm file: wrong offset after component 2"); if (read_data(binary_data, efficiency_factors, ByteOrder::little_endian) != Succeeded::yes) error("failed reading efficiency_factors from '%s'", full_data_file_name); + if (binary_data.tellg() != std::streampos(norm_parser.data_offset_each_dataset[3])) + error("Error reading ECAT8 norm file: wrong offset after component 3"); + if (read_data(binary_data, axial_effects, ByteOrder::little_endian) != Succeeded::yes) + error("failed reading axial_effects_factors from '%s'", full_data_file_name); if (scanner_ptr->get_type() == Scanner::Siemens_mMR) { @@ -468,14 +441,16 @@ MatrixFile* mptr = matrix_open(filename.c_str(), MAT_READ_ONLY, Norm3d); #endif -#if 1 - // to test pipe the obtained values into file - ofstream out_geom; + if (this->_write_components_to_file) + { + // to test pipe the obtained values into file + ofstream out_geom, out_axial; ofstream out_inter; ofstream out_eff; out_geom.open("geom_out.txt",ios::out); out_inter.open("inter_out.txt",ios::out); out_eff.open("eff_out.txt",ios::out); + out_axial.open("axial_out.txt",ios::out); for ( int i = geometric_factors.get_min_index(); i<=geometric_factors.get_max_index();i++) { @@ -505,7 +480,12 @@ MatrixFile* mptr = matrix_open(filename.c_str(), MAT_READ_ONLY, Norm3d); out_eff << std::endl<< std::endl; } -#endif + for ( int i = axial_effects.get_min_index(); i<=axial_effects.get_max_index();i++) + { + out_axial << axial_effects[i] << " " << std::endl; + } + + } #if 0 display(geometric_factors, "geo"); @@ -535,6 +515,13 @@ use_geometric_factors() const return this->_use_geometric_factors; } +bool +BinNormalisationFromECAT8:: +use_axial_effects_factors() const +{ + return this->_use_axial_effects_factors; +} + bool BinNormalisationFromECAT8:: use_crystal_interference_factors() const @@ -542,25 +529,9 @@ use_crystal_interference_factors() const return this->_use_crystal_interference_factors; } -#if 1 float BinNormalisationFromECAT8:: get_uncalibrated_bin_efficiency(const Bin& bin) const { - - // TODO disable when not HR+ or HR++ - /* - Additional correction for HR+ and HR++ - ====================================== - Modification of the normalisation based on segment number - Due to the difference in efficiency for the trues and scatter as the axial - angle increases - Scatter has a higher efficiency than trues when the axial angle is 0 (direct - planes) - As the axial angle increase the difference in efficiencies between trues and - scatter become closer - */ - const float geo_Z_corr = 1; - float total_efficiency = 0 ; @@ -599,7 +570,8 @@ get_uncalibrated_bin_efficiency(const Bin& bin) const { uncompressed_bin, detection_position_pair); - + const DetectionPosition<>& pos1 = detection_position_pair.pos1(); + const DetectionPosition<>& pos2 = detection_position_pair.pos2(); float lor_efficiency= 0.; /* @@ -664,12 +636,14 @@ get_uncalibrated_bin_efficiency(const Bin& bin) const { if (this->use_geometric_factors()) { lor_efficiency_this_pair *= -#ifdef SAME_AS_PETER - 1.F; -#else // this is 3dbkproj (at the moment) geometric_factors[geo_plane_num][uncompressed_bin.tangential_pos_num()]; -#endif } + if (this->use_axial_effects_factors()) + { + // Need to divide here + lor_efficiency_this_pair /= + find_axial_effects(pos1.axial_coord(), pos2.axial_coord()); + } lor_efficiency += lor_efficiency_this_pair; } @@ -682,31 +656,74 @@ get_uncalibrated_bin_efficiency(const Bin& bin) const { { view_efficiency += lor_efficiency; } + + total_efficiency += view_efficiency; } - - if (this->use_geometric_factors()) - { - /* z==bin.get_axial_pos_num() only when min_axial_pos_num()==0*/ - // for oblique plaanes use the single radial profile from segment 0 - -#ifdef SAME_AS_PETER - const int geo_plane_num = 0; - - total_efficiency += view_efficiency * - geometric_factors[geo_plane_num][uncompressed_bin.tangential_pos_num()] * - geo_Z_corr; -#else - total_efficiency += view_efficiency * geo_Z_corr; -#endif - } - else - { - total_efficiency += view_efficiency; - } } return total_efficiency; } -#endif + +void +BinNormalisationFromECAT8:: +construct_sino_lookup_table() +{ + const int num_rings = this->scanner_ptr->get_num_rings(); + this->sino_index=Array<2,int>(IndexRange2D(0, num_rings-1, + 0, num_rings-1)); + // fill with invalid index such that we can detect out-of-range + // see find_axial_effects + this->sino_index.fill(-1); + + auto proj_data_info_no_arccorr_ptr = + dynamic_cast(norm_proj_data_info_sptr.get()); + + if (!proj_data_info_no_arccorr_ptr) + error("BinNormalisationFromECAT8: internal error. Data should be of type ProjDataInfoCylindricalNoArcCorr"); + + this->num_Siemens_sinograms = proj_data_info_no_arccorr_ptr->get_num_non_tof_sinograms(); + + const auto segment_sequence = ecat::find_segment_sequence(*proj_data_info_no_arccorr_ptr); + Bin bin; + bin.tangential_pos_num()=0; + bin.view_num()=0; + std::vector > det_pos_pairs; + for (int Siemens_sino_index=0; Siemens_sino_index< this->num_Siemens_sinograms; ++Siemens_sino_index) + { + int z=Siemens_sino_index; + + for (std::size_t i=0; iget_num_axial_poss(bin.segment_num()); + if (z< num_ax_poss) + { + bin.axial_pos_num() = z; + proj_data_info_no_arccorr_ptr->get_all_det_pos_pairs_for_bin(det_pos_pairs, bin); + for (auto iter=det_pos_pairs.begin();iter!=det_pos_pairs.end(); ++iter) + { + sino_index[iter->pos1().axial_coord()][iter->pos2().axial_coord()] = Siemens_sino_index; + //sino_index[iter->pos2().axial_coord()][iter->pos1().axial_coord()] = Siemens_sino_index; + } + break; + } + else + { + z -= num_ax_poss; + } + } + } +} + +float +BinNormalisationFromECAT8:: +find_axial_effects(int ring1, int ring2) const +{ + const int Siemens_sino_index = sino_index[ring1][ring2]; + if (Siemens_sino_index<0) + return 1.F; + else + return axial_effects[Siemens_sino_index]; +} float diff --git a/src/recon_buildblock/BinNormalisationFromProjData.cxx b/src/recon_buildblock/BinNormalisationFromProjData.cxx index 060a59558..7af299da9 100644 --- a/src/recon_buildblock/BinNormalisationFromProjData.cxx +++ b/src/recon_buildblock/BinNormalisationFromProjData.cxx @@ -79,16 +79,25 @@ BinNormalisationFromProjData(const shared_ptr& norm_proj_data_ptr) Succeeded BinNormalisationFromProjData:: -set_up(const shared_ptr& exam_info_sptr, const shared_ptr& proj_data_info_ptr) +set_up(const shared_ptr& exam_info_sptr, const shared_ptr& proj_data_info_sptr) { - base_type::set_up(exam_info_sptr, proj_data_info_ptr); - - if (*(norm_proj_data_ptr->get_proj_data_info_sptr()) == *proj_data_info_ptr) + if (!base_type::set_up(exam_info_sptr, proj_data_info_sptr).succeeded()) + return Succeeded::no; + + const auto& norm_proj = *(norm_proj_data_ptr->get_proj_data_info_sptr()); + // complication: if the emission data is TOF but the norm is not, `apply()` et al. multiply all + // TOF bins with the non-TOF norm. So, we need to allow for that. + auto proj_to_check_sptr = proj_data_info_sptr; + if (!norm_proj.is_tof_data() && proj_data_info_sptr->is_tof_data()) + proj_to_check_sptr = proj_data_info_sptr->create_non_tof_clone(); + const auto& proj = *proj_to_check_sptr; + + if (norm_proj == proj) return Succeeded::yes; else { - const ProjDataInfo& norm_proj = *(norm_proj_data_ptr->get_proj_data_info_sptr()); - const ProjDataInfo& proj = *proj_data_info_ptr; + // Check if the emission data is "smaller" than the norm data (e.g. fewer segments) + // We will require equal tangential_pos ranges as `apply()` currently needs that. bool ok = (norm_proj >= proj) && (norm_proj.get_min_tangential_pos_num() ==proj.get_min_tangential_pos_num())&& @@ -98,7 +107,7 @@ set_up(const shared_ptr& exam_info_sptr, const shared_ptr& viewgrams, if (get_symmetries_used()->num_related_view_segment_numbers(basic_vs) != viewgrams.get_num_viewgrams()) - error("ForwardProjectByBin: forward_project called with incorrect related_viewgrams. Problem with symmetries!\n"); + error("ForwardProjectByBin: forward_project called with incorrect related_viewgrams (wrong number). Problem with symmetries!"); for (RelatedViewgrams::const_iterator iter = viewgrams.begin(); iter != viewgrams.end(); ++iter) { - ViewSegmentNumbers vs(iter->get_view_num(), iter->get_segment_num()); + ViewSegmentNumbers vs(iter->get_view_num(), iter->get_segment_num()); get_symmetries_used()->find_basic_view_segment_numbers(vs); + // TODOTOF find_basic_view_segment_numbers doesn't fill in timing_pos_num + vs.timing_pos_num() = basic_vs.timing_pos_num(); if (vs != basic_vs) error("ForwardProjectByBin: forward_project called with incorrect related_viewgrams. Problem with symmetries!\n"); } @@ -281,6 +283,8 @@ forward_project(RelatedViewgrams& viewgrams, { ViewSegmentNumbers vs(iter->get_view_num(), iter->get_segment_num()); get_symmetries_used()->find_basic_view_segment_numbers(vs); + // TODOTOF find_basic_view_segment_numbers doesn't fill in timing_pos_num + vs.timing_pos_num() = basic_vs.timing_pos_num(); if (vs != basic_vs) error("ForwardProjectByBin: forward_project called with incorrect related_viewgrams. Problem with symmetries!\n"); } diff --git a/src/recon_buildblock/GeneralisedObjectiveFunction.cxx b/src/recon_buildblock/GeneralisedObjectiveFunction.cxx index 1112d2a2d..ee2593908 100644 --- a/src/recon_buildblock/GeneralisedObjectiveFunction.cxx +++ b/src/recon_buildblock/GeneralisedObjectiveFunction.cxx @@ -41,7 +41,6 @@ set_defaults() this->prior_sptr.reset(); // note: cannot use set_num_subsets(1) here, as other parameters (such as projectors) are not set-up yet. this->num_subsets = 1; - use_tof = false; } template @@ -50,7 +49,6 @@ GeneralisedObjectiveFunction:: initialise_keymap() { this->parser.add_parsing_key("prior type", &prior_sptr); - this->parser.add_key("Use TOF information", &use_tof); } template @@ -82,11 +80,6 @@ set_up(shared_ptr const& target_data_ptr) return Succeeded::no; } - if (use_tof) - { - info("Time-Of-Flight reconstruction activated!"); - } - return Succeeded::yes; } @@ -183,14 +176,6 @@ get_num_subsets() const return this->num_subsets; } -template -bool -GeneralisedObjectiveFunction:: -get_tof_status() const -{ - return this->use_tof; -} - template double GeneralisedObjectiveFunction:: diff --git a/src/recon_buildblock/IterativeReconstruction.cxx b/src/recon_buildblock/IterativeReconstruction.cxx index 244542991..707461426 100644 --- a/src/recon_buildblock/IterativeReconstruction.cxx +++ b/src/recon_buildblock/IterativeReconstruction.cxx @@ -460,10 +460,10 @@ set_up(shared_ptr const& target_data_sptr) if (this->num_subsets<1 ) - { warning("number of subsets should be positive"); return Succeeded::no; } + { error("number of subsets should be positive"); return Succeeded::no; } if (is_null_ptr(this->objective_function_sptr)) - { warning("You have to specify an objective function"); return Succeeded::no; } + { error("You have to specify an objective function"); return Succeeded::no; } const int new_num_subsets = this->objective_function_sptr->set_num_subsets(this->num_subsets); @@ -474,24 +474,24 @@ set_up(shared_ptr const& target_data_sptr) this->num_subsets = new_num_subsets; } - if (this->objective_function_sptr->set_up(target_data_sptr) == Succeeded::no) - return Succeeded::no; - if (this->num_subiterations<1) - { warning("Range error in number of subiterations"); return Succeeded::no; } + { error("Range error in number of subiterations (has to be >= 1)"); return Succeeded::no; } if(this->start_subset_num<0 || this->start_subset_num>=this->num_subsets) - { warning("Range error in starting subset"); return Succeeded::no; } + { error("Range error in starting subset (has to be between 0 and num_subsets-1)"); return Succeeded::no; } if(this->save_interval<1 || this->save_interval>this->num_subiterations) - { warning("Range error in iteration save interval"); return Succeeded::no;} + { error("Range error in iteration save interval (has to be between 1 and num_subiterations-1)"); return Succeeded::no;} if (this->inter_iteration_filter_interval<0) - { warning("Range error in inter-iteration filter interval "); return Succeeded::no; } + { error("Range error in inter-iteration filter interval (has to be >= 0)"); return Succeeded::no; } if (this->start_subiteration_num<1) - { warning("Range error in starting subiteration number"); return Succeeded::no; } - + { error("Range error in starting subiteration number (has to be >= 1)"); return Succeeded::no; } + + if (this->objective_function_sptr->set_up(target_data_sptr) == Succeeded::no) + return Succeeded::no; + ////////////////// subset order // KT 05/07/2000 made randomise_subset_order int diff --git a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.cxx b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.cxx index b74620e85..94bd7c288 100644 --- a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.cxx +++ b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.cxx @@ -349,9 +349,9 @@ PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBinproj_data_info_sptr = this->list_mode_data_sptr->get_proj_data_info_sptr()->create_shared_clone(); +#if STIR_VERSION < 060000 if (this->get_max_segment_num_to_process() < 0) this->set_max_ring_difference(this->max_ring_difference_num_to_process); else @@ -472,6 +472,7 @@ PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBinnum_cache_files = 0; if(this->cache_lm_file) { info("Listmode reconstruction: Creating cache...", 2); diff --git a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx index a9ba5c88b..3ddc7cc81 100644 --- a/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx +++ b/src/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.cxx @@ -1,7 +1,7 @@ /* Copyright (C) 2000 PARAPET partners Copyright (C) 2000-2011, Hammersmith Imanet Ltd - Copyright (C) 2014, 2016-2022 University College London + Copyright (C) 2014, 2016-2023 University College London This file is part of STIR. SPDX-License-Identifier: Apache-2.0 AND License-ref-PARAPET-license @@ -550,13 +550,48 @@ sensitivity_uses_same_projector() const /*************************************************************** set_up() ***************************************************************/ +template +void +PoissonLogLikelihoodWithLinearModelForMeanAndProjData:: +ensure_norm_is_set_up(bool for_original_data) const +{ + + for_original_data = for_original_data || this->sensitivity_uses_same_projector(); + if (for_original_data) + { + if (!this->norm_already_setup || !this->latest_setup_norm_was_with_orig_data) + { + if (this->normalisation_sptr->set_up(proj_data_sptr->get_exam_info_sptr(), this->proj_data_sptr->get_proj_data_info_sptr()) == Succeeded::no) + error("Set_up of norm with original data failed."); + } + } + else + { + if (!this->norm_already_setup || this->latest_setup_norm_was_with_orig_data) + { + if (this->normalisation_sptr->set_up(proj_data_sptr->get_exam_info_sptr(), this->sens_proj_data_info_sptr) == Succeeded::no) + error("Set_up of norm with non-TOF data failed."); + } + } + this->norm_already_setup = true; + this->latest_setup_norm_was_with_orig_data = for_original_data; +} + +template +void +PoissonLogLikelihoodWithLinearModelForMeanAndProjData:: +ensure_norm_is_set_up_for_sensitivity() const +{ + this->ensure_norm_is_set_up(false); +} + template Succeeded PoissonLogLikelihoodWithLinearModelForMeanAndProjData:: set_up_before_sensitivity(shared_ptr const& target_sptr) { if (is_null_ptr(this->proj_data_sptr)) - error("you need to set the input data before calling set_up"); + error("you need to set the input data before calling set_up"); if (this->max_segment_num_to_process==-1) this->max_segment_num_to_process = @@ -565,12 +600,12 @@ set_up_before_sensitivity(shared_ptr const& target_sptr) if (this->max_segment_num_to_process > this->proj_data_sptr->get_max_segment_num()) { error("max_segment_num_to_process (%d) is too large", - this->max_segment_num_to_process); + this->max_segment_num_to_process); return Succeeded::no; } - this->max_timing_pos_num_to_process = - this->proj_data_sptr->get_max_tof_pos_num(); + this->max_timing_pos_num_to_process = + this->proj_data_sptr->get_max_tof_pos_num(); shared_ptr proj_data_info_sptr(this->proj_data_sptr->get_proj_data_info_sptr()->clone()); @@ -602,13 +637,15 @@ set_up_before_sensitivity(shared_ptr const& target_sptr) this->projector_pair_ptr->get_back_projector_sptr()->get_symmetries_used()->clone()); if (is_null_ptr(this->normalisation_sptr)) - { - error("Invalid normalisation object"); - return Succeeded::no; - } + { + error("Invalid normalisation object"); + return Succeeded::no; + } // we postpone calling setup_distributable_computation until we know which projectors we will use this->distributable_computation_already_setup = false; + // similar for norm + this->norm_already_setup = false; if (this->get_recompute_sensitivity()) { @@ -616,13 +653,13 @@ set_up_before_sensitivity(shared_ptr const& target_sptr) { this->sens_backprojector_sptr = projector_pair_ptr->get_back_projector_sptr(); this->sens_symmetries_sptr = this->symmetries_sptr; - if (this->normalisation_sptr->set_up(proj_data_sptr->get_exam_info_sptr(), proj_data_info_sptr) == Succeeded::no) - return Succeeded::no; + this->sens_proj_data_info_sptr = proj_data_info_sptr; } else { // sets non-tof backprojector for sensitivity calculation (clone of the back_projector + set projdatainfo to non-tof) auto pdi_non_tof_sptr = proj_data_info_sptr->create_non_tof_clone(); + this->sens_proj_data_info_sptr = pdi_non_tof_sptr; this->sens_backprojector_sptr.reset(projector_pair_ptr->get_back_projector_sptr()->clone()); if (auto sens_bp_pm_sptr = std::dynamic_pointer_cast(this->sens_backprojector_sptr)) { @@ -633,8 +670,6 @@ set_up_before_sensitivity(shared_ptr const& target_sptr) } this->sens_backprojector_sptr->set_up(pdi_non_tof_sptr, target_sptr); this->sens_symmetries_sptr.reset(this->sens_backprojector_sptr->get_symmetries_used()->clone()); - if (this->normalisation_sptr->set_up(proj_data_sptr->get_exam_info_sptr(), pdi_non_tof_sptr) == Succeeded::no) - return Succeeded::no; } } @@ -683,7 +718,8 @@ actual_compute_subset_gradient_without_penalty(TargetT& gradient, } if (!this->distributable_computation_already_setup) error("PoissonLogLikelihoodWithLinearModelForMeanAndProjData internal error: setup_distributable_computation not called (gradient calculation)"); - + if (add_sensitivity) + this->ensure_norm_is_set_up(); distributable_compute_gradient(this->projector_pair_ptr->get_forward_projector_sptr(), this->projector_pair_ptr->get_back_projector_sptr(), this->symmetries_sptr, @@ -725,6 +761,8 @@ actual_compute_objective_function_without_penalty(const TargetT& current_estimat } if (!this->distributable_computation_already_setup) error("PoissonLogLikelihoodWithLinearModelForMeanAndProjData internal error: setup_distributable_computation not called (function calculation)"); + this->ensure_norm_is_set_up(); + double accum=0.; distributable_accumulate_loglikelihood(this->projector_pair_ptr->get_forward_projector_sptr(), @@ -800,14 +838,9 @@ add_subset_sensitivity(TargetT& sensitivity, const int subset_num) const const int min_segment_num = -this->max_segment_num_to_process; const int max_segment_num = this->max_segment_num_to_process; -#if 1 shared_ptr sensitivity_this_subset_sptr(sensitivity.clone()); - shared_ptr sens_proj_data_sptr; // have to create a ProjData object filled with 1 here because otherwise zero_seg0_endplanes will not be effective - if (!this->sensitivity_uses_same_projector()) - sens_proj_data_sptr.reset(new ProjDataInMemory(this->proj_data_sptr->get_exam_info_sptr(), this->proj_data_sptr->get_proj_data_info_sptr()->create_non_tof_clone())); - else - sens_proj_data_sptr.reset(new ProjDataInMemory(this->proj_data_sptr->get_exam_info_sptr(), this->proj_data_sptr->get_proj_data_info_sptr())); + auto sens_proj_data_sptr = std::make_shared(this->proj_data_sptr->get_exam_info_sptr(), this->sens_proj_data_info_sptr); sens_proj_data_sptr->fill(1.0F); if (this->sensitivity_uses_same_projector() && (!this->distributable_computation_already_setup || !this->latest_setup_distributable_computation_was_with_orig_projectors)) @@ -836,10 +869,11 @@ add_subset_sensitivity(TargetT& sensitivity, const int subset_num) const this->distributable_computation_already_setup = true; this->latest_setup_distributable_computation_was_with_orig_projectors = false; } - if (!this->distributable_computation_already_setup) error("PoissonLogLikelihoodWithLinearModelForMeanAndProjData internal error: setup_distributable_computation not called (sensitivity calculation)"); + this->ensure_norm_is_set_up_for_sensitivity(); + distributable_sensitivity_computation(this->projector_pair_ptr->get_forward_projector_sptr(), this->sens_backprojector_sptr, this->sens_symmetries_sptr, @@ -863,66 +897,8 @@ add_subset_sensitivity(TargetT& sensitivity, const int subset_num) const std::transform(sensitivity.begin_all(), sensitivity.end_all(), sensitivity_this_subset_sptr->begin_all(), sensitivity.begin_all(), std::plus()); -#else - - // warning: has to be same as subset scheme used as in distributable_computation - for (int segment_num = min_segment_num; segment_num <= max_segment_num; ++segment_num) - { - //CPUTimer timer; - //timer.start(); - - for (int view = this->proj_data_sptr->get_min_view_num() + subset_num; - view <= this->proj_data_sptr->get_max_view_num(); - view += this->num_subsets) - { - const ViewSegmentNumbers view_segment_num(view, segment_num); - - if (!symmetries_sptr->is_basic(view_segment_num)) - continue; - this->add_view_seg_to_sensitivity(sensitivity, view_segment_num); - } - // cerr< -void -PoissonLogLikelihoodWithLinearModelForMeanAndProjData:: -add_view_seg_to_sensitivity(TargetT& sensitivity, const ViewSegmentNumbers& view_seg_nums) const -{ - int min_timing_pos_num = use_tofsens ? -this->max_timing_pos_num_to_process : 0; - int max_timing_pos_num = use_tofsens ? this->max_timing_pos_num_to_process : 0; - for (int timing_pos_num = min_timing_pos_num; timing_pos_num <= max_timing_pos_num; ++timing_pos_num) - { - RelatedViewgrams viewgrams = - this->proj_data_sptr->get_empty_related_viewgrams(view_seg_nums, - this->symmetries_sptr, false, timing_pos_num); - viewgrams.fill(1.F); - // find efficiencies - { - const double start_frame = this->frame_defs.get_start_time(this->frame_num); - const double end_frame = this->frame_defs.get_end_time(this->frame_num); - this->normalisation_sptr->undo(viewgrams, start_frame, end_frame); - } - // backproject - { - const int range_to_zero = - view_seg_nums.segment_num() == 0 && this->zero_seg0_end_planes - ? 1 : 0; - const int min_ax_pos_num = - viewgrams.get_min_axial_pos_num() + range_to_zero; - const int max_ax_pos_num = - viewgrams.get_max_axial_pos_num() - range_to_zero; - - this->sens_backprojector_sptr->back_project(sensitivity, viewgrams, min_ax_pos_num, max_ax_pos_num); - } - } - -} -#endif - template std::unique_ptr PoissonLogLikelihoodWithLinearModelForMeanAndProjData:: @@ -962,14 +938,11 @@ actual_add_multiplication_with_approximate_sub_Hessian_without_penalty(TargetT& } } + this->ensure_norm_is_set_up(); + shared_ptr symmetries_sptr( this->get_projector_pair().get_symmetries_used()->clone()); - const double start_time = - this->get_time_frame_definitions().get_start_time(this->get_time_frame_num()); - const double end_time = - this->get_time_frame_definitions().get_end_time(this->get_time_frame_num()); - this->get_projector_pair().get_forward_projector_sptr()->set_input(input); this->get_projector_pair().get_back_projector_sptr()->start_accumulating_in_new_target(); diff --git a/src/recon_buildblock/ProjMatrixByBin.cxx b/src/recon_buildblock/ProjMatrixByBin.cxx index ec4af0854..e402e9292 100644 --- a/src/recon_buildblock/ProjMatrixByBin.cxx +++ b/src/recon_buildblock/ProjMatrixByBin.cxx @@ -75,7 +75,6 @@ enable_tof(const shared_ptr& _proj_data_info_sptr, const boo if (v) { tof_enabled = true; - proj_data_info_sptr = _proj_data_info_sptr; gauss_sigma_in_mm = ProjDataInfo::tof_delta_time_to_mm(proj_data_info_sptr->get_scanner_ptr()->get_timing_resolution()) / 2.355f; r_sqrt2_gauss_sigma = 1.0f/ (gauss_sigma_in_mm * static_cast(sqrt(2.0))); } @@ -130,10 +129,13 @@ reserve_num_elements_in_cache(const std::size_t num_elems) void ProjMatrixByBin:: set_up( - const shared_ptr& proj_data_info_sptr, - const shared_ptr >& /*density_info_ptr*/ // TODO should be Info only + const shared_ptr& proj_data_info_sptr_v, + const shared_ptr >& density_info_sptr_v // TODO should be Info only ) { + this->proj_data_info_sptr = proj_data_info_sptr_v; + this->image_info_sptr.reset( + dynamic_cast* > (density_info_sptr_v->clone() )); if (is_cache_enabled()) { const int max_abs_tangential_pos_num = @@ -159,7 +161,6 @@ set_up( else { tof_enabled = false; - this->proj_data_info_sptr = proj_data_info_sptr; } this->cache_collection.recycle(); diff --git a/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx b/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx index 84d11dcba..01c587038 100644 --- a/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx +++ b/src/recon_buildblock/ProjMatrixByBinSPECTUB.cxx @@ -1,7 +1,7 @@ /* Copyright (C) 2013, Institute for Bioengineering of Catalonia Copyright (C) Biomedical Image Group (GIB), Universitat de Barcelona, Barcelona, Spain. - Copyright (C) 2013-2014, 2019, 2020 University College London + Copyright (C) 2013-2014, 2019, 2020, 2023 University College London Copyright (C) 2023 National Physical Laboratory This file is part of STIR. @@ -352,6 +352,15 @@ set_up( const VectorWithOffset radius_all_views = proj_Data_Info_Cylindrical->get_ring_radii_for_all_views(); + { + const auto max_radius = *std::max_element(radius_all_views.begin(), radius_all_views.end()); + const auto max_im_radius = std::max(vol.Xcmd2, vol.Ycmd2)*10; + if (max_im_radius > max_radius) + { + warning("Image radius (" + std::to_string(max_im_radius) + " is larger than max detector radius (" + + std::to_string(max_radius) + "). Are you sure this is correct? (Proceeding anyway)"); + } + } Rrad = new float [ wmh.prj.Nang ]; for ( int i = 0 ; i < wmh.prj.Nang ; i++ ) { // note: convert to cm for UB SPECT library @@ -374,7 +383,11 @@ set_up( //....PSF and collimator parameters ........................................ boost::algorithm::to_lower(psf_type); - if ( psf_type == "geometrical" ) wmh.do_psf = false; + if ( psf_type == "geometrical" ) + { + wmh.do_psf = false; + wmh.do_psf_3d = false; + } else{ wmh.do_psf = true; diff --git a/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx b/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx index 73a13481e..d305ee362 100644 --- a/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx +++ b/src/recon_buildblock/ProjMatrixByBinUsingRayTracing.cxx @@ -3,6 +3,7 @@ Copyright (C) 2000-2011, Hammersmith Imanet Ltd Copyright (C) 2013-2014, University College London Copyright (C) 2016, University of Hull + Copyright (C) 2013-2014, 2018, 2019, 2021, 2023 University College London Copyright 2017 ETH Zurich, Institute of Particle Physics and Astrophysics This file is part of STIR. @@ -271,18 +272,33 @@ set_up( const shared_ptr >& density_info_sptr_v // TODO should be Info only ) { - ProjMatrixByBin::set_up(proj_data_info_sptr_v, density_info_sptr_v); + auto image_info_ptr = dynamic_cast*> (density_info_sptr_v.get()); - image_info_sptr.reset( - dynamic_cast* > (density_info_sptr_v->clone() )); -// const VoxelsOnCartesianGrid * image_info_ptr = -// dynamic_cast*> (density_info_ptr.get()); + if (!image_info_ptr) + error("ProjMatrixByBinUsingRayTracing initialised with wrong type of DiscretisedDensity."); - if(is_null_ptr(image_info_sptr)) - error("ProjMatrixByBinUsingRayTracing initialised with a wrong type of DiscretisedDensity\n"); - - voxel_size = image_info_sptr->get_voxel_size(); - origin = image_info_sptr->get_origin(); + if (this->already_setup) + { + if (*this->proj_data_info_sptr == *proj_data_info_sptr_v && + this->voxel_size == image_info_ptr->get_voxel_size() && + this->origin == image_info_ptr->get_origin()) + { + CartesianCoordinate3D new_min_index; + CartesianCoordinate3D new_max_index; + image_info_ptr->get_regular_range(new_min_index, new_max_index); + if (this->max_index == new_max_index && + this->min_index == new_min_index) + { + info("ProjMatrixByBinUsingRayTracing::set_up skipped as already set-up with same characteristics.", 3); + } + return; + } + } + + ProjMatrixByBin::set_up(proj_data_info_sptr_v, density_info_sptr_v); + + voxel_size = image_info_ptr->get_voxel_size(); + origin = image_info_ptr->get_origin(); if (abs(origin.x())>.05F || abs(origin.y())>.05F) error("ProjMatrixByBinUsingRayTracing sadly doesn't support shifted x/y origin yet"); image_info_sptr->get_regular_range(min_index, max_index); @@ -300,14 +316,11 @@ set_up( if(sampling_distance_of_adjacent_LORs_xy/num_tangential_LORs > voxel_size.x() + 1.E-3 || sampling_distance_of_adjacent_LORs_xy/num_tangential_LORs > voxel_size.y() + 1.E-3) - warning("WARNING: ProjMatrixByBinUsingRayTracing used for pixel size (in x,y) " - "that is smaller than the bin size divided by num_tangential_LORs.\n" - "This matrix will completely miss some voxels for some (or all) views.\n"); - if(sampling_distance_of_adjacent_LORs_xy < voxel_size.x() - 1.E-3 || - sampling_distance_of_adjacent_LORs_xy < voxel_size.y() - 1.E-3) - warning("WARNING: ProjMatrixByBinUsingRayTracing used for pixel size (in x,y) " - "that is larger than the bin size.\n" - "Backprojecting with this matrix might have artefacts at views 0 and 90 degrees.\n"); + warning(boost::format("ProjMatrixByBinUsingRayTracing used for pixel size (x,y)=(%g,%g) " + "that is smaller than the central bin size (%g) divided by num_tangential_LORs (%d).\n" + "This matrix will completely miss some voxels for some (or all) views. It is therefore to best to increase " + "'number of rays in tangential direction to trace for each bin'.") + % voxel_size.x() % voxel_size.y() % sampling_distance_of_adjacent_LORs_xy % num_tangential_LORs); if (use_actual_detector_boundaries) { @@ -318,7 +331,7 @@ set_up( if (proj_data_info_cyl_ptr== 0) { warning("ProjMatrixByBinUsingRayTracing: use_actual_detector_boundaries" - " is reset to false as the projection data should be non-arccorected.\n"); + " is reset to false as the projection data should be non-arccorrected."); use_actual_detector_boundaries = false; } else @@ -348,7 +361,7 @@ set_up( if (proj_data_info_blk_ptr== 0) { warning("ProjMatrixByBinUsingRayTracing: use_actual_detector_boundaries" - " is reset to false as the projection data should be non-arccorected.\n"); + " is reset to false as the projection data should be non-arccorrected."); use_actual_detector_boundaries = false; } else @@ -365,7 +378,7 @@ set_up( if (!nocompression) { warning("ProjMatrixByBinUsingRayTracing: use_actual_detector_boundaries" - " is reset to false as the projection data as either mashed or uses axial compression\n"); + " is reset to false as the projection data as either mashed or uses axial compression."); use_actual_detector_boundaries = false; } } @@ -378,7 +391,7 @@ set_up( if (proj_data_info_blk_ptr== 0) { warning("ProjMatrixByBinUsingRayTracing: use_actual_detector_boundaries" - " is reset to false as the projection data should be non-arccorected.\n"); + " is reset to false as the projection data should be non-arccorrected."); use_actual_detector_boundaries = false; } else @@ -395,14 +408,14 @@ set_up( if (!nocompression) { warning("ProjMatrixByBinUsingRayTracing: use_actual_detector_boundaries" - " is reset to false as the projection data as either mashed or uses axial compression\n"); + " is reset to false as the projection data as either mashed or uses axial compression."); use_actual_detector_boundaries = false; } } } if (use_actual_detector_boundaries) - warning("ProjMatrixByBinUsingRayTracing: use_actual_detector_boundaries==true\n"); + info("ProjMatrixByBinUsingRayTracing: use_actual_detector_boundaries==true.", 3); } diff --git a/src/scatter_buildblock/CMakeLists.txt b/src/scatter_buildblock/CMakeLists.txt index abc5527cf..1a22da605 100644 --- a/src/scatter_buildblock/CMakeLists.txt +++ b/src/scatter_buildblock/CMakeLists.txt @@ -22,4 +22,4 @@ set(${dir_LIB_SOURCES} include(stir_lib_target) -target_link_libraries(scatter_buildblock recon_buildblock) +target_link_libraries(scatter_buildblock PUBLIC recon_buildblock) diff --git a/src/spatial_transformation_buildblock/CMakeLists.txt b/src/spatial_transformation_buildblock/CMakeLists.txt index a26196e43..0434d5836 100644 --- a/src/spatial_transformation_buildblock/CMakeLists.txt +++ b/src/spatial_transformation_buildblock/CMakeLists.txt @@ -23,4 +23,4 @@ set(${dir_LIB_SOURCES} include(stir_lib_target) -target_link_libraries(${dir} buildblock numerics_buildblock ) +target_link_libraries(${dir} PUBLIC buildblock numerics_buildblock ) diff --git a/src/swig/CMakeLists.txt b/src/swig/CMakeLists.txt index ecacf4900..6ebccab49 100644 --- a/src/swig/CMakeLists.txt +++ b/src/swig/CMakeLists.txt @@ -25,7 +25,7 @@ if(BUILD_SWIG_PYTHON OR BUILD_SWIG_OCTAVE OR BUILD_SWIG_MATLAB) FIND_PACKAGE(SWIG 3.0 REQUIRED) INCLUDE("${SWIG_USE_FILE}") - SET(CMAKE_SWIG_FLAGS -DSTART_NAMESPACE_STIR=\"namespace stir {\" -DEND_NAMESPACE_STIR=\"}\" -DSTIR_DEPRECATED="" -DSTIR_TOF=1) + SET(CMAKE_SWIG_FLAGS -DSTART_NAMESPACE_STIR=\"namespace stir {\" -DEND_NAMESPACE_STIR=\"}\" -DSTIR_DEPRECATED="" -DSTIR_TOF=1 -DSTIR_VERSION=${VERSION}) # Manually add some flags (really should get this from elsewhere) if (HAVE_LLN_MATRIX) SET(CMAKE_SWIG_FLAGS ${CMAKE_SWIG_FLAGS} -DHAVE_LLN_MATRIX) @@ -156,8 +156,8 @@ if(BUILD_SWIG_PYTHON) endif() endif() SWIG_WORKAROUND(${SWIG_MODULE_stir_REAL_NAME}) - SWIG_LINK_LIBRARIES(stir ${STIR_LIBRARIES} ${STIR_Python_dependency}) - target_link_libraries(${SWIG_MODULE_stir_REAL_NAME} ${OpenMP_EXE_LINKER_FLAGS}) + SWIG_LINK_LIBRARIES(stir PUBLIC ${STIR_LIBRARIES} ${STIR_Python_dependency}) + target_link_libraries(${SWIG_MODULE_stir_REAL_NAME} PUBLIC ${OpenMP_EXE_LINKER_FLAGS}) CONFIGURE_FILE(./pyfragments.swg ./ COPYONLY) set(PYTHON_DEST ${CMAKE_INSTALL_PREFIX}/python CACHE PATH "Destination for python module") @@ -207,7 +207,7 @@ if (BUILD_SWIG_OCTAVE) endif() SET_TARGET_PROPERTIES(${SWIG_MODULE_stiroct_REAL_NAME} PROPERTIES SUFFIX ${OCTAVE_SUFFIX} PREFIX "${OCTAVE_PREFIX}") SWIG_WORKAROUND(${SWIG_MODULE_stiroct_REAL_NAME}) - SWIG_LINK_LIBRARIES(stiroct ${STIR_LIBRARIES} ${OCTAVE_LIBRARIES}) + SWIG_LINK_LIBRARIES(stiroct PUBLIC ${STIR_LIBRARIES} ${OCTAVE_LIBRARIES}) # add OCTAVE_INCFLAGS to swig-generated file only, not to all files as # 1) we don't need it at the moment 2) we'd need to change from -Ibla to bla @@ -242,9 +242,9 @@ if (BUILD_SWIG_MATLAB) LINK_FLAGS "${Matlab_CXXLINKER_FLAGS}" FOLDER "Matlab") SWIG_WORKAROUND(${SWIG_MODULE_stirMATLAB_REAL_NAME}) - target_link_libraries(${SWIG_MODULE_stirMATLAB_REAL_NAME} ${OpenMP_EXE_LINKER_FLAGS}) + target_link_libraries(${SWIG_MODULE_stirMATLAB_REAL_NAME} PUBLIC ${OpenMP_EXE_LINKER_FLAGS}) - SWIG_LINK_LIBRARIES(stirMATLAB ${STIR_LIBRARIES} ${Matlab_LIBRARIES}) + SWIG_LINK_LIBRARIES(stirMATLAB PUBLIC ${STIR_LIBRARIES} ${Matlab_LIBRARIES}) include_directories(${Matlab_INCLUDE_DIRS}) # disabled, as currently set via add_definitions in main CMakeLists.txt diff --git a/src/swig/stir.i b/src/swig/stir.i index c4cd86ec7..ea5f1f18c 100644 --- a/src/swig/stir.i +++ b/src/swig/stir.i @@ -1,6 +1,6 @@ /* Copyright (C) 2011-07-01 - 2012, Kris Thielemans - Copyright (C) 2013, 2014, 2015, 2018 - 2022 University College London + Copyright (C) 2013, 2014, 2015, 2018 - 2023 University College London Copyright (C) 2022 National Physical Laboratory Copyright (C) 2022 Positrigo This file is part of STIR. @@ -30,12 +30,14 @@ #include // for size_t #include #include + #ifdef SWIGOCTAVE // TODO terrible work-around to avoid conflict between stir::error and Octave error // they are in conflict with eachother because we need "using namespace stir" below (swig bug) #define __stir_error_H__ #endif +#include "stir/stream.h" // to get access to stream output for STIR types for ADD_REPR etc #include "stir/num_threads.h" #include "stir/find_STIR_config.h" @@ -61,6 +63,9 @@ #include "stir/copy_fill.h" #include "stir/ProjDataInterfile.h" + #include "stir/Radionuclide.h" + #include "stir/RadionuclideDB.h" + #include "stir/DataSymmetriesForViewSegmentNumbers.h" #include "stir/recon_buildblock/BinNormalisationFromProjData.h" #include "stir/recon_buildblock/BinNormalisationFromAttenuationImage.h" @@ -112,6 +117,8 @@ #include "stir/HUToMuImageProcessor.h" #include "stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.h" +#include "stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeData.h" +#include "stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.h" #include "stir/OSMAPOSL/OSMAPOSLReconstruction.h" #include "stir/OSSPS/OSSPSReconstruction.h" #include "stir/recon_buildblock/ForwardProjectorByBinUsingProjMatrixByBin.h" @@ -567,6 +574,10 @@ %ignore *::ask_parameters; %ignore *::create_shared_clone; %ignore *::read_from_stream; +%ignore *::get_data_ptr; +%ignore *::get_const_data_ptr; +%ignore *::release_data_ptr; +%ignore *::release_const_data_ptr; #if defined(SWIGPYTHON) %rename(__assign__) *::operator=; @@ -890,6 +901,7 @@ namespace std { // Macros for adding __repr_()_ for Python and disp() for MATLAB // example usage: ADD_REPR(stir::ImagingModality, %arg($self->get_name())); + // second argument piped to stream, so could be a std::string, but also another type %define ADD_REPR(classname, defrepr) %extend classname { @@ -900,18 +912,20 @@ namespace std { #if defined(SWIGPYTHON) std::string __repr__() { - std::string repr = ""; + return s.str(); } #endif #if defined(SWIGMATLAB) void disp() { - std::string repr = "<" + "classname::"; - repr += defrepr; - repr += ">"; + std::stringstream s; + s << ""; mexPrintf(repr.c_str()); } #endif @@ -1000,11 +1014,26 @@ namespace std { /* Parse the header files to generate wrappers */ //%include "stir/shared_ptr.h" %include "stir/Succeeded.h" +ADD_REPR(stir::Succeeded, %arg($self->succeeded() ? "yes" : "no")); + %include "stir/NumericType.h" %include "stir/ByteOrder.h" -%include "stir/DetectionPosition.h" %include "stir_coordinates.i" + +#if 0 + // TODO enable this in STIR version 6 (breaks backwards compatibility +%attributeref(stir::LORInAxialAndNoArcCorrSinogramCoordinates, float, z1); +%attributeref(stir::LORInAxialAndNoArcCorrSinogramCoordinates, float, z2); +%attributeref(stir::LORInAxialAndNoArcCorrSinogramCoordinates, float, beta); +%attributeref(stir::LORInAxialAndNoArcCorrSinogramCoordinates, float, phi); +#else +%ignore *::z1() const; +%ignore *::z2() const; +%ignore *::beta() const; +%ignore *::phi() const; +#endif +%ignore *::check_state; %include "stir/LORCoordinates.h" %template(FloatLOR) stir::LOR; @@ -1084,6 +1113,12 @@ namespace std { %shared_ptr(stir::FanProjData); %shared_ptr(stir::GeoData3D); +%ignore operator<<; +%ignore operator>>; +%ignore stir::DetPairData::operator()(const int a, const int b) const; +%ignore stir::DetPairData3D::operator()(const int a, const int b) const; +%ignore stir::FanProjData::operator()(const int, const int, const int, const int) const; +%ignore stir::GeoData3D::operator()(const int, const int, const int, const int) const; %include "stir/ML_norm.h" %shared_ptr(stir::InvertAxis); diff --git a/src/swig/stir_exam.i b/src/swig/stir_exam.i index 4a14d9c90..fb2e49102 100644 --- a/src/swig/stir_exam.i +++ b/src/swig/stir_exam.i @@ -16,6 +16,7 @@ %shared_ptr(stir::TimeFrameDefinitions); %shared_ptr(stir::ImagingModality); %shared_ptr(stir::PatientPosition); +%shared_ptr(stir::RadionuclideDB); %shared_ptr(stir::Radionuclide); %shared_ptr(stir::ExamInfo); %shared_ptr(stir::ExamData); @@ -24,6 +25,7 @@ %include "stir/ImagingModality.h" %include "stir/PatientPosition.h" %include "stir/Radionuclide.h" +%include "stir/RadionuclideDB.h" %include "stir/ExamInfo.h" %include "stir/ExamData.h" diff --git a/src/swig/stir_objectivefunctions.i b/src/swig/stir_objectivefunctions.i index 2f5ccd8b8..7876271f1 100644 --- a/src/swig/stir_objectivefunctions.i +++ b/src/swig/stir_objectivefunctions.i @@ -32,11 +32,16 @@ %shared_ptr(stir::GeneralisedObjectiveFunction); %shared_ptr(stir::PoissonLogLikelihoodWithLinearModelForMean); +%shared_ptr(stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData); %shared_ptr(stir::RegisteredParsingObject, stir::GeneralisedObjectiveFunction, stir::PoissonLogLikelihoodWithLinearModelForMean >); +%shared_ptr(stir::RegisteredParsingObject, + stir::GeneralisedObjectiveFunction, + stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData >); %shared_ptr(stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData); +%shared_ptr(stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin); %shared_ptr(stir::SqrtHessianRowSum); @@ -48,6 +53,8 @@ %include "stir/recon_buildblock/GeneralisedObjectiveFunction.h" %include "stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMean.h" %include "stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndProjData.h" +%include "stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeData.h" +%include "stir/recon_buildblock/PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin.h" %include "stir/recon_buildblock/SqrtHessianRowSum.h" @@ -57,6 +64,7 @@ %template (GeneralisedObjectiveFunction3DFloat) stir::GeneralisedObjectiveFunction; //%template () stir::GeneralisedObjectiveFunction; %template (PoissonLogLikelihoodWithLinearModelForMean3DFloat) stir::PoissonLogLikelihoodWithLinearModelForMean; +%template (PoissonLogLikelihoodWithLinearModelForMeanAndListModeData3DFloat) stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData; // TODO do we really need this name? // Without it we don't see the parsing functions in python... @@ -65,14 +73,19 @@ stir::GeneralisedObjectiveFunction, stir::PoissonLogLikelihoodWithLinearModelForMean >; +%template(RPPoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin3DFloat) stir::RegisteredParsingObject, + stir::GeneralisedObjectiveFunction, + stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeData >; + %template (PoissonLogLikelihoodWithLinearModelForMeanAndProjData3DFloat) stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData; +%template (PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin3DFloat) stir::PoissonLogLikelihoodWithLinearModelForMeanAndListModeDataWithProjMatrixByBin; %inline %{ template stir::PoissonLogLikelihoodWithLinearModelForMeanAndProjData * ToPoissonLogLikelihoodWithLinearModelForMeanAndProjData(stir::GeneralisedObjectiveFunction *b) { return dynamic_cast*>(b); -} + } %} %template(ToPoissonLogLikelihoodWithLinearModelForMeanAndProjData3DFloat) ToPoissonLogLikelihoodWithLinearModelForMeanAndProjData; diff --git a/src/swig/stir_projdata.i b/src/swig/stir_projdata.i index 78cfb548c..9911d1e39 100644 --- a/src/swig/stir_projdata.i +++ b/src/swig/stir_projdata.i @@ -1,6 +1,6 @@ /* Copyright (C) 2011-07-01 - 2012, Kris Thielemans - Copyright (C) 2013, 2014, 2015, 2018 - 2022 University College London + Copyright (C) 2013, 2014, 2015, 2018 - 2022, 2023 University College London Copyright (C) 2022 National Physical Laboratory This file is part of STIR. @@ -34,16 +34,47 @@ %shared_ptr(stir::Sinogram); %shared_ptr(stir::Viewgram); -%newobject stir::Scanner::get_scanner_from_name; -%include "stir/Scanner.h" +%shared_ptr(stir::DetectionPosition); +%shared_ptr(stir::DetectionPositionPair); + +%attributeref(stir::DetectionPosition, unsigned int, tangential_coord); +%attributeref(stir::DetectionPosition, unsigned int, axial_coord); +%attributeref(stir::DetectionPosition, unsigned int, radial_coord); +%include "stir/DetectionPosition.h" +#ifdef STIR_TOF +ADD_REPR(stir::DetectionPosition, %arg(*$self)) +#endif +%template(DetectionPosition) stir::DetectionPosition; + +%attributeref(stir::DetectionPositionPair, stir::DetectionPosition, pos1); +%attributeref(stir::DetectionPositionPair, stir::DetectionPosition, pos2); +%include "stir/DetectionPositionPair.h" +#ifdef STIR_TOF + //ADD_REPR(stir::DetectionPositionPair, %arg(*$self)) +#endif +%template(DetectionPositionPair) stir::DetectionPositionPair; -%attributeref(stir::Bin, int, segment_num); +%attributeref(stir::SegmentIndices, int, segment_num); +#ifdef STIR_TOF +%attributeref(stir::SegmentIndices, int, timing_pos_num); +#endif +%attributeref(stir::ViewgramIndices, int, view_num); +%attributeref(stir::SinogramIndices, int, axial_pos_num); %attributeref(stir::Bin, int, axial_pos_num); -%attributeref(stir::Bin, int, view_num); %attributeref(stir::Bin, int, tangential_pos_num); %attributeref(stir::Bin, int, timing_pos_num); +%attributeref(stir::Bin, int, time_frame_num); %attribute(stir::Bin, float, bin_value, get_bin_value, set_bin_value); +%include "stir/SegmentIndices.h" +%include "stir/ViewgramIndices.h" +%include "stir/SinogramIndices.h" %include "stir/Bin.h" +#ifdef STIR_TOF +ADD_REPR(stir::Bin, %arg(*$self)) +#endif + +%newobject stir::Scanner::get_scanner_from_name; +%include "stir/Scanner.h" %newobject stir::ProjDataInfo::ProjDataInfoGE; %newobject stir::ProjDataInfo::ProjDataInfoCTI; diff --git a/src/swig/stir_reconstruction.i b/src/swig/stir_reconstruction.i index efc992cbc..863c2ba0f 100644 --- a/src/swig/stir_reconstruction.i +++ b/src/swig/stir_reconstruction.i @@ -15,7 +15,7 @@ \author Kris Thielemans \author Markus Jehl */ -%rename (get_inter_iteration_filter) *::get_inter_iteration_filter_sptr; +%ignore *::get_inter_iteration_filter_sptr; %rename (get_subset_sensitivity) *::get_subset_sensitivity_sptr; %rename (set_objective_function) *::set_objective_function_sptr; %ignore *::get_objective_function_sptr; // we have it without _sptr in C++ diff --git a/src/swig/test/python/test_buildblock.py b/src/swig/test/python/test_buildblock.py index d615c0da8..2d2472b7d 100755 --- a/src/swig/test/python/test_buildblock.py +++ b/src/swig/test/python/test_buildblock.py @@ -4,7 +4,7 @@ # py.test test_buildblock.py -# Copyright (C) 2013 University College London +# Copyright (C) 2013, 2023 University College London # This file is part of STIR. # # SPDX-License-Identifier: Apache-2.0 @@ -187,14 +187,32 @@ def test_zoom_image(): assert abs(zoomed_image[ind]-1./(zoom))<.001 def test_Scanner(): - s=Scanner.get_scanner_from_name("ECAT 962") - assert s.get_num_rings()==32 - assert s.get_num_detectors_per_ring()==576 - #l=s.get_all_names() - #print s + scanner=Scanner.get_scanner_from_name("ECAT 962") + assert scanner.get_num_rings()==32 + assert scanner.get_num_detectors_per_ring()==576 + #l=scanner.get_all_names() + #print scanner # does not work #for a in l: # print a + scanner=Scanner.get_scanner_from_name("SAFIRDualRingPrototype") + scanner.set_scanner_geometry("BlocksOnCylindrical") + scanner.set_up() + d=DetectionPosition(1,1,0) + c=scanner.get_coordinate_for_det_pos(d) + d2=DetectionPosition(); + s=scanner.find_detection_position_given_cartesian_coordinate(d2, c) + assert s.succeeded() + assert d==d2 + +def test_Radionuclide(): + modality = ImagingModality(ImagingModality.PT) + db = RadionuclideDB() + r = db.get_radionuclide(modality, "^18^Fluorine") + assert abs(r.get_half_life() - 6584) < 1 + modality = ImagingModality(ImagingModality.NM) + r = db.get_radionuclide(modality, "^99m^Technetium") + assert abs(r.get_half_life() - 6.0058*3600) < 10 def test_Bin(): segment_num=1; @@ -214,6 +232,8 @@ def test_Bin(): assert abs(bin.bin_value-bin_value)<.01; bin=Bin(segment_num, view_num, axial_pos_num, tangential_pos_num, bin_value); assert abs(bin.bin_value-bin_value)<.01; + bin.time_frame_num=3; + assert bin.time_frame_num==3; def test_ProjDataInfo(): s=Scanner.get_scanner_from_name("ECAT 962") diff --git a/src/test/test_interpolate_projdata.cxx b/src/test/test_interpolate_projdata.cxx index 4589ad338..06e43fac9 100644 --- a/src/test/test_interpolate_projdata.cxx +++ b/src/test/test_interpolate_projdata.cxx @@ -58,6 +58,7 @@ class InterpolationTests : public RunTests void check_symmetry(const SegmentBySinogram& segment); void compare_segment(const SegmentBySinogram& segment1, const SegmentBySinogram& segment2, float maxDiff); + void compare_segment_shape(const SegmentBySinogram& shape_segment, const SegmentBySinogram& test_segment, int erosion); }; void InterpolationTests::check_symmetry(const SegmentBySinogram& segment) @@ -147,6 +148,46 @@ void InterpolationTests::compare_segment(const SegmentBySinogram& segment check_if_less(sumAbsDifference, maxDiff, "difference between segments is larger than expected"); } +void InterpolationTests::compare_segment_shape(const SegmentBySinogram& shape_segment, const SegmentBySinogram& test_segment, int erosion) +{ + auto maxTestValue = test_segment.find_max(); + // compute difference and compare against empirically found value from visually validated sinograms + auto sumVoxelsOutsideMask = 0U; + for (auto axial = test_segment.get_min_axial_pos_num(); axial <= test_segment.get_max_axial_pos_num(); axial++) + { + for (auto view = test_segment.get_min_view_num(); view <= test_segment.get_max_view_num(); view++) + { + for (auto tang = test_segment.get_min_tangential_pos_num(); tang <= test_segment.get_max_tangential_pos_num(); tang++) + { + if (test_segment[axial][view][tang] < 0.1 * maxTestValue) + continue; + + // now go through the erosion neighbourhood of the voxel to see if it is near a non-zero voxel + bool isNearNonZero = false; + for (auto axialShape = std::max(axial - erosion, test_segment.get_min_axial_pos_num()); + axialShape <= std::min(axial + erosion, test_segment.get_max_axial_pos_num()); axialShape++) + { + for (auto viewShape = std::max(view - erosion, test_segment.get_min_view_num()); + viewShape <= std::min(view + erosion, test_segment.get_max_view_num()); viewShape++) + { + for (auto tangShape = std::max(tang - erosion, test_segment.get_min_tangential_pos_num()); + tangShape <= std::min(tang + erosion, test_segment.get_max_tangential_pos_num()); tangShape++) + { + if (shape_segment[axialShape][viewShape][tangShape] > 0) + isNearNonZero = true; + } + } + } + if (isNearNonZero == false) + sumVoxelsOutsideMask++; + } + } + } + + // confirm that the difference is smaller than an empirically found value + check_if_equal(sumVoxelsOutsideMask, 0U, "there were non-zero voxels outside the masked area"); +} + void InterpolationTests::scatter_interpolation_test_blocks() { info("Performing symmetric interpolation test for BlocksOnCylindrical scanner"); @@ -172,8 +213,9 @@ void InterpolationTests::scatter_interpolation_test_blocks() // define a cylinder precisely in the middle of the FOV, such that symmetry can be used for validation auto emission_map = VoxelsOnCartesianGrid(*downsampled_proj_data_info, 1); - auto cylinder = EllipsoidalCylinder(40, 80, 80, CartesianCoordinate3D(emission_map.get_max_z()/2 * emission_map.get_grid_spacing()[1], 0, 0)); + auto cylinder = EllipsoidalCylinder(41, 80, 80, CartesianCoordinate3D(emission_map.get_max_z()/2 * emission_map.get_grid_spacing()[1], 0, 0)); cylinder.construct_volume(emission_map, CartesianCoordinate3D(1, 1, 1)); + write_to_file("downsampled_cylinder_map", emission_map); // project the cylinder onto the downsampled scanner proj data auto pm = ProjMatrixByBinUsingRayTracing(); @@ -316,7 +358,7 @@ void InterpolationTests::scatter_interpolation_test_blocks_asymmetric() interpolated_proj_data.write_to_file("interpolated_sino_asym.hs"); // compare to ground truth - compare_segment(interpolated_proj_data.get_segment_by_sinogram(0), full_size_model_sino.get_segment_by_sinogram(0), 4697); + 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_cyl_asymmetric() @@ -384,7 +426,7 @@ void InterpolationTests::scatter_interpolation_test_cyl_asymmetric() interpolated_proj_data.write_to_file("interpolated_sino_cyl_asym.hs"); // compare to ground truth - compare_segment(interpolated_proj_data.get_segment_by_sinogram(0), full_size_model_sino.get_segment_by_sinogram(0), 17500); + compare_segment_shape(full_size_model_sino.get_segment_by_sinogram(0), interpolated_proj_data.get_segment_by_sinogram(0), 2); } void diff --git a/src/utilities/UPENN/CMakeLists.txt b/src/utilities/UPENN/CMakeLists.txt index 54fe8d133..c0d6a0d77 100644 --- a/src/utilities/UPENN/CMakeLists.txt +++ b/src/utilities/UPENN/CMakeLists.txt @@ -12,5 +12,5 @@ set(${dir_EXE_SOURCES} include(stir_exe_targets) -target_link_libraries(conv_UPENN_projdata_to_STIR ${UPENN_libgeom} +target_link_libraries(conv_UPENN_projdata_to_STIR PUBLIC ${UPENN_libgeom} ${UPENN_liblor} ${UPENN_libimagio} ${UPENN_libimagio++} -lboost_system)