Summary of changes in STIR release 4.0

Overall summary

This version is 95% backwards compatible with STIR 3.0 for the user (see below). Developers might need to make code changes as detailed below. In addition, this release adds a lot of new functionality to STIR. Of course, there is also the usual code-cleanup and improvements to the documentation.

This release contains mainly code written by

• Kris Thielemans (UCL) for general overview, fixes, help and maintenance
• Nikos Efthimiou (UCL and Hull) for the scatter update, list mode reconstruction, ROOT support etc
• Richard Brown (UCL) for improving code structure and safety (e.g. projectors, normalisation, set_up checks), IO additions, handling of dynamic data
• Daniel Deidda (Leeds and NPL) for the PLS prior, the HKEM algorithm, and list mode and projection functionalities for SPECT
• Ashley Gillman (CSIRO AeHRC) for patient orientation and geometrical support

Various bits were tested by (in alphabetical order, all UCL except where indicated): Ottavia Bertolli, Matthias Ehrhardt, Pawel Markiewicz, Richard Manber, Alaleh Rashidnasab, Benjamin Thomas (CIRC ASTAR and UCL), Yu-jung Tsai, Vesna Cuplov.

Patch release info

• 4.0.0 released 22/04/2020
• 4.0.1 released 28/04/2020
  • PR 513 (suppress warnings with clang)
  • issue 515, PR 516 (compilation problem with JSON present and BUILD_EXECUTABLES=OFF)
  • issue 520 (siemens mMR unlisting segfault)
  • typos in UsersGuide
• 4.0.2 released 22/05/2020
  • PR 540 (KOSMAPOSL stability improvements and fixes in naming of class member function)

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

Changes breaking backwards compatibility from a user-perspective

• major changes in the scatter code, see the section on scatter updates
• ECAT6 is no longer supported
• minor changes when using Python/MATLAB (see below)
• the need for a new z zoom keyword for SPECT reconstructions (see the updated sample .par files)
• changed image orientation when reading/writing images via ITK. LPS coordinate info is now preserved.
• replaced Gaussian filters in experimental with better, but incompatible, versions in buildblock. See below for more information.
• Compared to previous versions of STIR, the SPECT reconstructed images are now inverted along the x axis to match DICOM and nifti (and PET) orientation.
• InterfileDynamicDiscretisedDensityOutputFileFormat changed start/end key from Interfile Output File Format to Interfile Output File Format Parameters to match other Interfile types.
• Renamed InterfileParametricDensityOutputFileFormat to InterfileParametricDiscretisedDensityOutputFileFormat and the corresponding keys.
• >fwdtest and bcktest are now only installed when BUILD_TESTING is ON. Use forward_project and backward_project instead.

Bug fixes

• bug fixed in writing of SPECT interfile headers:
  • images were always written as PET. This was mostly harmless, except when using the image then as input to forward_project, in which case the output projection data was written with a PET Interfile header, resulting in incorrect geometry etc etc.
  • for projection data, "direction of rotation" was set incorrectly. This affected output of forward_project (and stir_math -s). Reconstruction was not affected though.
• fix conversion from raw GATE using conv_GATE_raw_ECAT_projdata_to_interfile(it was using the wrong segment-order and data-type).
• the function value of the quadratic prior value was inconsistent with its gradient. It is now divided by 2. This bug did not affect OSMAPOSL nor OSSPS, but did affect the value of the objective function being printed out. This bug would cause major problems if you developed your own reconstruction algorithms that used the objective function value. (bug found by Yu-Jung Tsai).
• minor bug fix in FBP3DRP. The determination of the rmin/rmax of the missing data was incorrect for negative delta (the range was too narrow). In principle, this could have caused problems, but only when a projector was used where the segment symmetry is switched off.
• When using the inNormalisationFromAttenuationImage constructor that takes an image, the code would rescale the image of the caller.
  This bug did not affect command line utilities, but did affect python/matlab or use of the C++ library.

New functionality

• Major changes in the (PET) scatter code, see the section on scatter updates
• Extra reconstruction capabilities:
  • KOSMAPOSL STIR now allows to use the hybrid kernel method described in the following paper:
    Deidda D, et al., Hybrid PET-MR list-mode kernelized expectation maximization reconstruction,
    Inverse Problems, 2019, DOI: 10.1088/1361-6420/ab013f.
    KOSMAPOSL also allows to use multiple anatomical images. The method is described in
    Deidda D, et al., Multiplexing Kernelized Expectation Maximization Reconstruction for PET-MR
    NSS/MIC proceedings, 2018, DOI: 10.1109/NSSMIC.2018.8824312. KOSMAPOSL allows to use openMP for the kernel calculation. KOSMAPOSL allows reconstruction without any anatomical image
  • Image-based PSF-modelling in PET can be achieved by using projectors that first smooth data image before forward projection or after back projection. However, this is currently somewhat difficult to set-up as it needs specifying the filter twice. Note that the SPECT projector does allow sinogram-based PSF modelling.
  • Added the anatomical parallel level sets (PLS) prior. If an uniform anatomical image is provided it will work as a (smoothed) total variation (TV) prior. See
    Ehrhardt MJ, et al., PET Reconstruction With an Anatomical MRI Prior Using Parallel Level Sets
    IEEE Trans Med Imaging. 2016 Sep;35(9):2189-2199. doi: 10.1109/TMI.2016.2549601. Epub 2016 Apr 14.
  • Reconstruction of listmode data has been corrected. LM-OSEM can be performed but only geometrical subsets are currently supported.
  • Reconstruction parameter files now can use a new optional keyword z zoom to allow changing the spacing from the default. It is currently only recommended for SPECT (where it should be set to 0.5 due current work-arounds on how to read SPECT data into STIR). PET projectors are likely to fail when using the non-default z zoom.
  • Indirect parametric reconstruction (Patlak) no longer requires you to pass a pre-existing template file. If one is supplied, it will be used; else the parametric image will be created based on the dynamic image.
  • Normalisation for SPECT has been added. Can be used with lm_to_projdata and during reconstruction. Currently provides the option to apply uniformity and decay correction.
• Parallelisation:
  • Sensitivity calculation parallelized using distributed worker mechanism.
  • This release considerably expands OpenMP support for multi-threaded execution. Previously, only FBP2D was multi-threaded and it underperformed. Now, any reconstruction code using distributable_computation() (including OSMAPOSL and OSSPS) and the scatter estimation code is multi-threaded. See
    Kris Thielemans, Vesna Cuplov, Benjamin A. Thomas, Multi-threaded image reconstruction of 3D PET sinogram data with STIR,
    proc. Fully 3D Conference, Newport, RI, USA, Jun 2015.
    By default, we use about 90% of your cores, but you can change that by setting the OMP_NUM_THREADS environment variable before you start your STIR work (see the wiki). Please note that you cannot enable both OPENMP and MPI.
• Filters:
  • Added SeparableGaussianArrayFilter and SeparableGaussianImageFilter (based on old code in experimental which has now been removed). This filter constructs the kernel by sampling a Gaussian. It is therefore different from the Metz filter with power 0, as that is constructed by sampling a Gaussian in frequency space (and can therefore have negative values in image space).
• IO:
  • Meta-data:
    • we now keep track of time frames and patient orientation for images, in particular making sure that reconstructions initialise this information from the projection data (if it was present there). In the next version, providing this information will likely become a requirement.
    • It is now highly recommended to specify imaging modality and patient orientation in your Interfile headers. If you don't, STIR will assume you know what you're doing and write a warning (but we might not cover all cases). For patient orientation, we will assume HFS.
  • The orientation of images read/written via ITK has changed. It should now be correct if the patient was in HFS position. Unfortunately, for most file formats (such as NifTI) this can't be verified due to lack of metadata in the format. If the files are in DICOM format, as long as the file contains the correct metadata, other patient positions are supported (currently HFS, FFS, HFP or FFP).
  • Siemens file formats:
    • Support the parsing of SIEMENS Interfile files natively, except for images. You still have to use convertSiemensInterfileToSTIR.sh for those.
    • List-mode classes for Siemens ECAT8 (32-bit), tested on the mMR only though.
    • Normalisation class for Siemens ECAT8, allowing for the virtual, tested on the mMR only though. This needed 2 extra functions in Scanner to know how many virtual crystals there are.
  • GATE compatibility layer. ROOT files from cylindricalPET and ECAT simulated scanners can be imported as listmode files. Some simple options are provided, like exclude scattered or randoms. A new header file (recommended extension .hroot is necessary for ROOT files. Check the examples.
  • Ability to read dynamic and parametric interfiles with InterfileDynamicDiscretisedDensityInputFileFormat and InterfileParametricDiscretisedDensityInputFileFormat, respectively.
  • The creation of MultipleDataSetHeader. This enables the user to read a set of filenames from a text file, referred to as a multi file. See examples/samples/multi_file.txt for an example.
  • List-Mode functionalities are now extended to SPECT. However, we currently cannot distribute any software for SPECT scanners due to confidentiality clauses.
  • CListModeData got a new function get_total_number_of_events(), BUT currently only CListModeDataFromROOT supports it.
  • When writing ECAT7 files, patient orientation and time frame info were not set.
  • Support for Energy Resolution and energy windows in scanners and projection data has been added. It is only required when running the scatter simulation or estimation, and should then be set in the input file. For example, for a STIR-Interfile header, you should use the following keywords:

        number of energy windows:=1
        energy window lower level[1]:=425
        energy window upper level[1]:=650
        Scanner parameters:=
          Scanner type := ...
          Energy resolution := 0.145
          Reference energy (in keV) := 511
          ...
        End scanner parameters:=
          

    Only 1 energy window is currently supported (if there's more than one, only the first will be read).
    Warning: Unfortunately, GitHub versions between version 3.0 and 4.0 used

            energy window lower level:=425
            energy window upper level:=650
          

    This will still be read ok, as long as these keywords are preceded by

            !version of keys := STIR3.0
          

    The new (vectored) version is compatible with the Interfile standard.
• Python (and MATLAB) interface
  • when BUILD_DOCUMENTATION is ON, we now use doxy2swig to pass the doxygen documentation to Python/MATLAB.
  • added GeneralisedPoissonNoiseGenerator
  • added the Shape classes
  • added the GeneralisedPrior classes
  see also the developer section below
• New or changed utilities:
  • list_image_info now has options to list exam info, or only part of the information (as list_projdata_info).
  • list_image_values now has options to print coordinates in LPS as opposed to indices, print in CSV, or suppress the title row.
  • list_lm_countrates is a new list mode utility to output total counts per specified time interval.
  • Extract single images from dynamic image with executable extract_dynamic_images.
  • list_ROI_values now has the option to save the image name , max and min as new columns, and they will appear in the same row as the roi name, mean value and so on.
  • New functionality added to zoom_image. The user can now specify different types of zoom: preserve_values, preserve_projections, preserve_sum.
  • generate_image now has optional keywords for modality, patient position and timing.
  • create_projdata_template now allows even span. In particular, span 2 corresponds to how GE stores data. See the STIR glossary for more information.
  • ctac_to_mu_values is a new utility that converts a CT image to linear attenuation coefficients. This application depends on the nlohmann JSON library.
  • manip_image now has an extra option that allows to substitute nan with zero.
  • New utility invert_axis to invert x, y or z axis of an image.
  • lm_to_projdata_with_random_rejections is a new addition that makes possible to randomly reduce the number of event in a list mode file. The events are binned into a sinogram.
  • Triple energy Window (TEW) for (SPECT) scatter correction is now available. Needs availability of different sinograms acquired at different energy windows.
  • Functions and executables for ML estimation of singles/randoms/normalisation factors for PET, written mostly by Kris Thielemans, then at Hammersmith Imanet Ltd. This code is kindly contributed by GE Healthcare.. These are only minimally documented at the present stage unfortunately.
    • Hogg, K. Thielemans, T. Spinks, N. Spyrou,
      Maximum-Likelihood Estimation of Normalisation Factors for PET,
      proc. of IEEE Medical Imaging Conf. 2001, vol. 4 pp. 2065 - 2069. DOI: 10.1109/NSSMIC.2001.1009231
    • D. Hogg, K. Thielemans, S. Mustafovic and T.J. Spinks, A Study of Bias for various Iterative Reconstruction Methods in PET,
      Proc. of IEEE Medical Imaging Conf. 2002, vol. 3, pp. 1519 - 1523. DOI: 10.1109/NSSMIC.2002.1239610
• add Verbosity level. It can be set from interactive languages, and in reconstructions with the verbosity keyword in the parameter file.

Changed functionality

• Projection- and normalisation-related classes have had changes to make them much safer to use in an interactive environment such as Python. The methods will now check if set_up() has been called and if the images/projection data have the same characteristics as those used when the object was set-up. (Actually, projection data can be "smaller", e.g. less segments).
• create_projdata_template now returns an extra segment pair if the maximum ring difference does not coincide with a "full" segment. This change was to accomodate how GE uses even span.
• The num_events_to_store variable in LmToProjData has been upgraded to long int allow processing of more events in list mode files.
• calculate_attenuation_coefficients now defaults to using the ray-tracing matrix as projector. Results should be identical though. This means it will now automatically handle cases where there are less symmetries in the data, e.g. for scanners with odd number of views.
  It has a new switch --NOPMRT which reverts to the previous behaviour (i.e. using the old ray tracer on-the-fly projector).
• upsample_and_fit_single_scatter previously aborted if there were not enough counts (after weighting with the masking used for the tail-fit) in the scatter estimate for one sinogram. Now we just write a warning and set it to 1. For consistency, we now use scale factor 1 if all sinograms had zero counts (after weighting), while we were setting it 0 before. As this should never occur, and the scale factor should be irrelevant then, this should not affect any results.
  (Actual change is in get_scale_factors_per_sinogram).
• Changed the return value of DynamicProjData::read_from_file from a raw ptr to a unique_ptr.
• Output template for get_dynamic_images_from_parametric_images is now optional. Writing of the output is no longer restricted to ecat7.
• Code that was hosted in src/local and include/local/stir is now in src/experimental and include/stir_experimental, respectively. An extra CMake parameter has also been added, STIR_ENABLE_EXPERIMENTAL, which is OFF by default. Enabling this parameter will include the experimental code in the build.
• The DICOM series reading restrictions in read_file_itk have been relaxed. The requirement for the same acquisition, trigger or frame times in a series is no longer enforced. This allows 3-D DICOM volumes to be read correctly but may cause issues for dynamic PET or CT data.
• Data processors (such as motion or smoothing) can be used prior to forward projections and following back projections. This supercedes functionality provided by PresmoothingForwardProjectorByBin and PostsmoothingBackProjectorByBin, which can now be considered deprecated.

Python (and MATLAB) interface

• renamed get_scanner_ptr() and get_proj_data_info_ptr() in the get_scanner() and get_proj_data_info()
• expose Bin member functions as variables (or attributes). As opposed to bin.view_num(), write bin.view_num. You can now also do bin.view_num=1.

Build system

The build system now only uses CMake. We recommend to use CMake version 3.1 or later on Linux/MacOSX and 3.10.2 or later on Windows, but we try to accomodate older versions of CMake.

• MATLAB support now uses CMake's own FindMatlab module, as opposed to our own code, which was out-of-date. Unfortunately, this uses different CMake variables. In particular, MATLAB_ROOT is now ignored. Use Matlab_ROOT_DIR instead. In addition, CMake launches MATLAB to find out some settings, which on some systems can be very slow.
• Also find the Python executable. This allows specifying your python version by defining PYTHON_EXECUTABLE when running CMake.
• added (advanced) options STIR_USE_BOOST_SHARED_PTR and STIR_NO_UNIQUE_PTR. These default to Off unless Boost believes that your compiler does not have a working std::shared_ptr or std::unique_ptr (and in particular if you didn't enable C++11 or higher). You should normally not have to change these options.
• provide CMake option to disable ITK usage even if it's found.
• fixes related to OpenMP
• Changes to CMake build files, allowing others to import STIR via find_package. To allow this, CMake now creates a config.h file which sets preprocessor variables etc (as opposed to passing these on the compiler command line).
• added BUILD_TESTING and BUILD_EXECUTABLES options to reduce build time if you do not need those files
• install include files
• if doxygen is found, added the target RUN_DOXYGEN for building the Doxygen documentation (thanks to Jannis Fischer). This target is always built if BUILD_DOCUMENTATION is ON.
• Added a BUILD_TESTS build target so that make BUILD_TESTS (or whatever is appropriate for your build system) will build all tests (not the other executables or Python etc libraries) .

Known problems

See the Github issues for more information.

• SPECT related
  • Do NOT use current PET projectors on SPECT Data. You will have to choose the "SPECT UB" projector in all cases (see the sample .par files).
  • The quantitative scale of reconstructed SPECT images is generally different than for PET reconstructions. It probably also behaves differently than the PET projectors when using different voxel sizes.
  • Images are still written as PET Interfile (but this is ignored when reading).
  • At present, the collimator model varies the width of the PSF linearly with distance. This was experimentally verified with some collimators (for normal distances to the collimator). However, we need to add the possibility to use different models.
  • list_projdata_info still reports SPECT data as if it's a PET scanner, but you can see all necessary geometric info from it anyway.
  • create_projdata_template only works for PET scanners. For SPECT, you could take the sample Interfile header, modify it to your needs (and create an empty file for the binary data pointed to by the Interfile header).
• Building with shared libraries is currently broken on many (all?) systems due to inter-dependencies between libraries.
• This release is known to fail when both MPI and shared libraries are enabled. You will see linking errors in all non_MPI executables.
• The incremental interpolating back-projector still gives problems in the centre of the images on some architectures/compiler/optimisation settings. We recommend to use the ray-tracing matrix (with additional rays) instead.
• When using GRAPHICS=X (i.e. on Unix/Linux/MacOSX) the window displaying the bitmaps is sometimes empty. This might depend on your X server settings. This affects manip_image, manip_projdata and display_projdata. We have currently no work-around (aside from using GRAPHICS=PGM or stir_write_pgm).

Minor bug fixes

• fix HighResWallClockTimer on Linux and Win32 which caused reporting wrong timings in certain situations
• LmToProjData crashed when a template was used which was incompatible with the unlisting (e.g. a arc-corrected template for most cylindrical scanners).
• listmode files could not be re-read when the end of the file was reached. This could have affected cases where the number of segments was limited and all events were read while unlisting the file. (e.g. in lm_to_projdata)

Documentation changes

• Added documentation on new features
• The Wiki expands more and more, so use it in addition to the provided PDFs.

recon_test_pack changes

• updated version number and added some clarification to the README.txt
• Additional test scripts, e.g. for GATE produced ROOT files, listmode reconstruction and scatter simulation, have been added.

Other changes to tests

• added tests for the ITK IO routines and Interfile reading of dynamic/parametric images/projection data

Scatter updates

• estimate_scatter is now renamed to simulate_scatter, as that is what it does. Its parameter format changed as well. See an example in examples/PET_simulation. In particular:
  • start and end keywords are different
  • energy window and resolution information is now (only) taken from the template projection data.
  • density_image_filename keywords are renamed to attenuation_image_filename.
  • template_proj_data_filename is renamed to template_projdata_filename for consistency with other parameter files.
  Functionality and results are identical to the previous version.
• Somewhat confusingly, estimate_scatter is a new executable that actually estimates scatter from measured data. It performs the iterative process first suggested by Watson et al. (similar to Ollinger et al.). An example of this as in examples/Siemens-mMR/scatter.
• In the github version after STIR 3.0, we included a shell script for easier scatter estimation. That script still exists for backwards compatibility, but is essentially obsolete. You should use estimate_scatter instead.

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

Major bugs fixed

• see above

Backward incompatibities

• The base class for list mode has been generalised from CListMode to ListMode. SPECTListMode for SPECT and CListMode for PET are derived classes. All the classes and utilities for list mode now use the ListMode class (ListModeData, ListRecord, ListEvent etc. accordingly) and new developments using list mode need to use the base class as well.
• The ScatterEstimation class has been renamed to ScatterSimulation. ScatterEstimation is a new class that performs an iterative process to estimate scatter.
• STIR now uses std::unique_ptr as opposed to std::auto_ptr, as recommended for C++-11 and later (unless the CMake option STIR_NO_UNIQUE_PTR is set). This affects calls to read_from_file and DataSymmetriesForBins::find_symmetry_operation_from_basic_bin etc. You should normally just be able to do a simple find-and-replace.
  Note that if std::unique_ptr is not defined for your compiler (and compilation flags), STIR will try to work-around this (see src/include/stir/unique_ptr.h).
• Make stir::shared_ptr to be std::shared_ptr (stick to boost::shared_ptr if the CMake option STIR_USE_BOOST_SHARED_PTR is set, i.e. normally when BOOST_NO_CXX11_SMART_PTR is defined). We now also import dynamic_pointer_cast and static_pointer_cast into the stir namespace such that you can write code that doesn't need to know which shared_ptr is used. Finally, use the preprocessor define MAKE_SHARED as opposed to std::make_shared.
• The image must be set prior to using the forward and back projectors. For the forward projector this is done with set_input. For the back projector, this is done during the set_up, but thenstart_accumulating_in_new_target must be called each time the user wishes to start back projecting into a zero-filled image. The result of the back projection is then returned via get_output. Backwards compatibility can be maintained with the CMake switch, STIR_PROJECTORS_AS_V3. The user should not notice these changes unless they have developed their own projectors. These changes were made so that using data processors (such as motion or smoothing) prior to the forwards projection or following the back projection will be more efficient.
• Projection- and BinNormalisation-related derived classes now have to call set_up() of their base class. Otherwise, you will get run-time errors saying that set_up() was not called.
• ProjDataInfo::ProjDataInfoCTI (superseded by ProjDataInfo::construct_proj_data_info) now returns an extra segment pair if the maximum ring difference does not coincide with a "full" segment. This change was to accomodate how GE uses even span.
• Classes derived from CListModeData now need to call set_proj_data_info_sptr and can then use get_proj_data_info_sptr(). This replaces the scanner_sptr protected member (which no longer exists). This will break all derived classes.
• ProjMatrixByBin now uses get_symmetries_sptr() (replaces get_symmetries_ptr()).
• The inheritance between Object, ParsingObject and RegisteredObject have changed. Specifically:
  • Object inherits from ParsingObject.
  • Object no longer needs the method parameter_info().
  • Classes that inherit from RegisteredObject no longer need to also inherit from ParsingObject.
  • The class AddParser has been removed. The default third argument of RegisteredParsingObject has been changed accordingly.
  • Object has been renamed as RegisteredObjectBase. This is because the class only exists for the parser to get information about the class.
• KeyParser has changed to be safer for parsing of vectorised keys. However, the "simple" add_key functions now create non-vectorised keys. Use add_vectorised_key otherwise.

New functionality

• write_to_file is a new addition that makes it easier to write using the default output file format.
• FilePath is a new addition that provides a thin cross-platform layer for file and folder operations. In addition to the features already found in utilities it can create new paths and perform checks on paths about accessibility and writting permissions.
• Added apply/undo methods to BinNormalisation to process whole ProjData
• extra functionality for ProjData:
  • added fill() functions (with float or with another ProjData)
  • added get_proj_data_info_sptr() (also to Sinogram et al)
  • write_to_file moved from ProjDataInMemory to ProjData
• Information on a single energy window has been added in the ExamInfo.
• ExamData is a new base class for ProjData, CListModeData, DiscretisedDensity and the dynamic/gated variations.
• Moved functions from the poisson_noise utility to a new class GeneralisedPoissonNoiseGenerator.
• Added the ability to KeyParser to parse long int numbers. Also add remove_key.
• Introduce operator>= for ProjDataInfo.
• New functionality added to zoom_image. The user can now specify different types of zoom: preserve_values, preserve_projections, preserve_sum.

Other code changes

• We removed all using std::vector etc statements from .h,.inl,.txx files. This means you might have to modify your own code accordingly.
• Changed ProjMatrixByBin to have separate caches per view-segment number to allow much larger axial/tangential_pos range and for finer-grained OPENMP locking.
  WARNING: breaks backwards compatibility as ProjMatrixByBin::set_up now has to be called by set_up in each derived class.
• info() and warning() really should no longer be used in the printf style, i.e. warning("Some message %d", x) . In the next version, they will only accept string-type argument (and optional verbosity). Use boost::format.
• get_proj_data_info_ptr() and get_scanner_ptr() should no longer be used as they can create problems with the lifetime of objects. Use get..._sptr instead.
• ClistEvent::is_valid_template is a new (pure) virtual function which is called by LmToProjData to check if the list mode format can handle the template. All derived classes of CListEvent will need to implement it.
• We removed SeparableMetzArrayFilter::discrete_fourier_transform(). Now stir::fourier is used. No need to comply with the Licence of "Numerical Recipies" anymore.