v5.0.0

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

Overall summary

• At least C++-11 is now required. We are not aware of any problems with more recent versions of C++.
• Initial support for PET scanners with block detectors or even generic location of detectors (less tested feature), taking into account appropriate geometry. This is currently kept independent of the cylindrical scanner modelling used normally, but this will be changed in a future version. See PR #577. This work is described in
  P. Khateri, J. Fischer, W. Lustermann, C. Tsoumpas, and G. Dissertori,
  Implementation of cylindrical PET scanners with block detector geometry in STIR,
  EJNMMI Physics, vol. 6, no. 1, p. 15, Jul. 2019, doi: 10.1186/s40658-019-0248-9.
  V. Dao, E. Mikhaylova, M. L. Ahnen, J. Fischer , K. Thielemans, C. Tsoumpas,
  Image Reconstruction for PET Scanners With Polygonal Prism Geometry Using STIR, proc. IEEE MIC 2021
  This code was initially developed by Parisa Khateri and Michael Roethlisberger (ETH), and further improved and tested by Viet Dao (Leeds), Daniel Deidda (NPL) and Kris Thielemans (UCL and ASC) with funding by the University of Leeds and Positrigo AG.
• View Offset Support enabled for the PET scanners, contributed by Palak Wadhwa (Leeds Univ), allowing fixing rotation of reconstructed images, see PR 181.
• Maximum Likelihood estimation of normalisation factors now includes estimation of geometric factors in 3D, see PR 619. This code was mostly contributed by Tahereh Niknejad (work completed at University of Lisbao, Portugal and PETsys Electronics, together with Kris Thielemans, UCL), as well as capabilities to work with virtual crystals, see PR 833 and PR 949. See the following proceedings for the initial work
  Tahereh Niknejad, Stefaan Tavernier, Joao Varela, and Kris Thielemans, Validation of 3D Model-Based Maximum-Likelihood Estimation of Normalisation Factors for Partial Ring Positron Emission Tomography In 2016 IEEE Nuclear Science Symposium, Medical Imaging Conference and Room-Temperature Semiconductor Detector Workshop (NSS/MIC/RTSD), 1–5. DOI: 10.1109/NSSMIC.2016.8069577.
• A new ProjDataInfoSubsetByView class has been implemented which allows for subsets of projection data. Previously, subset projections could only be stored in fully sampled ProjData. This is still possible, the implementation is backward compatible, and subset reconstructions such as OSMAPSOL have not been changed. However, there is now an option to create subsets with ProjData::get_subset which only store a subset of views and so are more memory efficient when subsets might be handled separately. ProjDataInfoSubsetByView currently does not support file I/O.
• ROOT files produced by GATE can now be interpreted using "virtual crystals", see PR #617.
• A calibration factor can be applied in the normalisation with a new BinNormalisation derived class BinNormalisationWithCalibration; Interfile now reads/writes isotope name and calibration factor. See PR 672. This code was contributed by Daniel Deidda (NPL) and Kris Thielemans (UCL)
• Radionuclide information is saved into a JSON configuration file for a number of popular radionuclide for PET and SPECT. This can make operations like decay correction easier and accurate. To be able to use the database we need to provide the isotope name used in the acquisition. For SPECT this is extracted from the Dicom header and saved in the interfile. For PET ways to read this information need to be implemented for each scanner/vendor. The new classes are new Radionuclide and RadionuclideDB; This code was contributed by Daniel Deidda (NPL) and Kris Thielemans (UCL)
• KOSMAPOSL (HKEM) allows now to freeze the iterative part of the kernel at a specific subiteration. Code contributed by Daniel Deidda (NPL), Kris Thielemans (UCL) and Ashley Gillman (CSIRO)
• It is now possibly to call the Georg's Schramm parallelproj projectors (using either OpenMP or CUDA), see PR 817, contributed by Kris Thielemans (UCL) with Georg Schramm (KUL).
• OpenMP loop scheduling changed to use dynamic instead of runtime, resulting in faster performance.

Of course, there is also the usual code-cleanup and improvements to the documentation.

This release contains mainly code written by Kris Thielemans (UCL), Richard Brown (UCL), Parisa Khateri and Michael Roethlisberger (ETHZ), Robert Twyman (UCL), Daniel Deidda (NPL), Tahereh Nikjenad (PETsys), Palak Wadhwa (Univ of Leeds), Viet Ahn Dao (Univ of Leeds), Ashley Gillman (CSIRO), Georg Schramm (KUL), Markus Jehl (Positrigo), Gefei Chen (Univ of Macao)

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

Changes breaking backwards compatibility from a user-perspective

• View Offset Support enabled for the PET scanners.

  • For the scanners that have non-zero intrinsic azimuthal tilt angle, reconstructed images will now get rotated.
  • If you use view mashing, in previous STIR versions images were rotated according to half the number of mashed views. This is now corrected.

  WARNING This means that reconstructed PET images will not be identical when either the scanner has non-zero view-offset, or view-mashing is on.
  To reflect this change, Scanner::get_default_intrinsic_tilt() has been renamed to get_intrinsic_azimuthal_tilt().
  Note: start angle was already supported for SPECT
  Backward compatibility for reconstructed images can be achieved by setting the STIR_LEGACY_IGNORE_VIEW_OFFSET CMake option to ON. (However, copied sinogram data will then always have the offset set to 0, in contrast to earlier versions of STIR).

• PoissonLogLikelihood Hessian methods have been corrected to reflect the concavity of the function. Hessian vector products now return non-positive voxels, if the vector (input) is non-negative. STIR usages of these methods (OSSPS and SqrtHessianRowSum) have been updated to see no effective change in functionality. However, calling the Hessian vector product methods, via python etc., will result in a change in functionality as the sign of the output voxel values has changed.

• Python/MATLAB wrappers no longer have ProjDataInfo::ProjDataInfoCTI, use ProjDataInfo::construct_proj_data_info instead. However, this now returns an object of the appropriate (derived) type as opposed to just ProjDataInfo. This should be transparent, except apparently for testing equality of objects.

Bug fixes

• There was an inconsistency between log-likelihood function and its gradient when use_subset_sensitivities was false and num_subsets greater than 1. Fixing this isssue means that images reconstructed with OSSPS are different from previous versions of STIR when the above conditions are met. See Issue #873 and associated PR #893.
• Parametric image reconstruction with POSMAPOSL could lead to zeroes being introduced gradually during reconstruction. See Issue #906 and associated PR #978.

New functionality

• Capability to model block and generic PET detectors was added. This is currently limited to span=1, no view mashing and the ray-tracing matrix (single LOR). It is enabled by specifying appropriate keywords in the Interfile header of the projection data.

        Scanner parameters:=
        normal parameters...
        Scanner geometry (BlocksOnCylindrical/Cylindrical/Generic)  := BlocksOnCylindrical
        Distance between crystals in axial direction (cm)           := 0.22
        Distance between crystals in transaxial direction (cm)      := 0.22
        Distance between blocks in axial direction (cm)             := 0.22
        Distance between blocks in transaxial direction (cm)        := 3.3
        end scanner parameters:=
  

  Scatter and normalisation code are still pending changes.

• Georg's Schramm parallelproj is an implementation of the Joseph algorithm for projection. If it has been installed in your system, and you tell CMake where to find it (`parallelproj_DIR=/wherever/lib/cmake`), the STIR user is now able to select an additional projector, called Parallelproj. This will use the CUDA version if available, otherwise will fall-back to the OpenMP version. Check the new sample files in examples/samples and the section in the User's Guide.

• The (still preliminary) code for Maximum Likelihood estimation of normalisation factors now includes estimation of geometric factors in 3D as well

• The (also preliminary) ProjDataInfoSubsetByView is backwards compatible and allows a new, memory efficient method for subset projections. Subset projections are created with with ProjData::get_subset, which is the only addition to the ProjData class (additions are mostly in the aforementioned new Info class). ProjDataInfoSubsetByView can be used to set up projectors, which will project only the subset. ProjDataInfoSubsetByView currently does not support file I/O. Therefore, the interface is not yet accessible through parameter files, only through C++ interface.
  PR #969

• calculate_attenuation_coefficients utility now accepts an optional forward projection parameter file in the same format as the forward_project utility. Example usage has been added to src/recon_test_pack/simulate_data.sh.

• ROOT files produced by GATE can now be interpreted using ``virtual crystals", i.e. by inserting ``dummy" crystals before converting to cylindrical geometry (as is done on many Siemens scanners). LmToProjData and list mode reconstructions will then put the LORs more accurately (at least when the ``virtual crystals" roughly correspond to gaps between blocks). See the update .hroot file in examples/samples.
  Warning: if you use the originating system to specify the scanner, this will be automatically enabled. (If you do not want this, set it to User_defined_scanner and specify all values).
  PR #617

• Moved most code from the ctac_to_mu_values.cxx utility to a new HUToMuImageProcessor class (derived from DataProcessor. This allows combining it with filtering etc, also from in Python. It means that the ctac_to_mu_values.cxx utility itself is now obsolete, but it still exists in this version. (Note that this functionality depends on an external JSON library.)

• Addition of a new utility create_multi_header which can be used to create a single header pointing to several files (e.g. one image per time frame). The header uses the (STIR-specific) multi format.

• Addition of a new utility extract_single_images_from_parametric_image to get each parametric image in a single file.

• generate_image parameter files now support the originating system keyword.

• list_image_info now works for dynamic images, with a new --per-volume option to list min/max/sum for every volume.

• SSRB now has the option of taking a template sinogram.

• KOSMAPOSL (HKEM) allows now to freeze the iterative part of the kernel at a specific subiteration. This can be set in the parameter file through the keyword freeze_iterative_kernel_at_subiter_num

• BinNormalisationWithCalibration is a new class derived from BinNormalisation. This class allows to apply calibration factor from the scanner and save the information into the interfile header. To use this, a specific BinNormalisation class should be derived from BinNormalisationWithCalibration, the information about the calibration factor read, and the function get_uncalibrated_bin_efficiencies() needs to be overwritten. Note that also the isotope name and branching ratio can be set here (the latter will need to be set according to the isotope by reading into a radionuclide database, see below). BinNormalisationSPECT already reads Calibration factor and isotope name and apply the calibration factor read from interfile. Since not all SPECT scanner do quantitative reconstruction an option of setting the calibration factor from the parameter file is added. Factors for BinNormalisationECA8/ECAT7 and GEHDF5 are set to one at the moment as we need to double check on the meaning of cross calibration factor for ECAT and find out how to read them for GEHDF5. Documentation for this is pending and will be added in the following PR.

• Radionuclide is a new class which contains radionuclide information such as halflife, branching ratio, energy etc. This class allows to propagate infrmation trough the reconstruction as it is now a member in the ExamInfo. It is populated from the information extracted from RadionuclideDB which allows to create a Radionuclide object from the information stored in the radionuclide_info.json in the STIR configuration directory. Different vendors, as we saw with SPECT, may have different standards for the isotope name. Three different ones were observed, therefore a lookup table has been added to allow the use of the data base for any different isotope name format. The lookup table is a JSON file stored in radionuclide_names.json.

• Logcosh Prior.

• SAFIR input file format now supports a Gaussian LOR randomization, which is applied on the endpoints of each LOR when sorting the listmode events into virtual scanner projection data.

• Scanner has now a new static member get_names_of_predefined_scanners returning a list of names.

• Exposed compute_total_ROI_values and ROIValues to Python/MATLAB via swig. Added python example ROI_demo.py to demonstrate usage.

• Additional C++ demo demonstrating the use the objective function and gradient ascent optimisation, see examples/src/demo4_obj_fun.cxx.

• A few functions related to TOF were added to ProjDataInfo and Scanner for future compatibility. However, this release still does not support TOF. (Check the GitHub site for the relevant PR).

• Exposed BinNormalisation, ListModeData, LmToProjData, ProjectorByBinPair, ScatterSimulation, ScatterEstimation and a couple of support classes and functions to Python/MATLAB via swig. Added python example listmode_demo.py to demonstrate usage.

• More setter functions on ScatterEstimation.

Changed functionality

• OpenMP loop scheduling changed to use dynamic instead of runtime. On some machines, this was causing slower projection operations due to a static scheduler being selected by default. See Issue #935 for more details.
• Many operations with ProjDataInMemory are now much faster (it now uses an underlying 1D array).
• ParametricDiscretisedDensity objects can now have an ExamInfo object with multiple time frames (corresponding to the time frames of the data where the parametric image is derived from). In some cases, there could only be a single time frame (start to end of the study).
• find_ML_normfactors3D and apply_ML_normfactors3D can be used for scanner with virtual crystal (Only verified on mCT). Contributed by Gefei Chen, see PR 833.
• If there are multiple buckets specified in the interfile header, we increase the symmetry size to a bucket in the find_ML_normfactors3D. Otherwise, we use a block. The use of an argument --for-symmetry-per-block in the find_ML_normfactors3D utility will use the symmetry size of a block.
• apply_patlak_to_images no longer uses an existing file as a template for the dynamic image but will overwrite it.
• ROOT file I/O improvements. An entire entry's tree is no longer loaded by default and instead individual branches are loaded as needed. ROOT file event processing is now up to 4x faster. In addition, there is now an option to check energy window information (defaulting to on). Futhermore, added functionality to exclude true and unscattered event types from list mode processing.

Build system

• We now require CMake at least version 3.1, although we highly recommended to use a very recent version of CMake to avoid problems with libraries or compilers which are more recent than your CMake version.
• At least C++-11 is now required. We are not aware of any problems with most recent versions of C++. When building, you can change the C++ version by setting CMAKE_CXX_STANDARD, see the CMake documentation for supported values.
  When importing STIR's STIRConfig.cmake via find_package(STIR), your compiler will be set to use at least C++-11 (via CMake's target_compile_features).
  Note that some external libraries that STIR depends on (such as ROOT) might increase the required C++ version, depending on how they were built.
• CERN's ROOT library is now preferentially found by searching for its own exported ROOTConfig.cmake. Set the CMake variable ROOT_DIR accordingly. Older behaviour relying on ROOTSYS and root-config will be deprecated in a future version.
• If you have CMake 3.14 or more recent, Python is found using the FindPython module, as opposed to the deprecated FindPythonLibs. FindPythonInterp modules. You can set the Python_EXECUTABLE variable to a particular version (see the CMake doc for other hints you can provide).( If the old PYTHON_EXECUTABLE is set, we use it to initialise Python_EXECUTABLE).

Known problems

Minor bug fixes

• Removed boolean state from SAFIR ListMode factory, to enable reading more than one listmode file in the same runtime environment.

Documentation changes

• Added documentation on new features
• Also check the wiki in addition to the provided PDFs.

recon_test_pack changes

• updated version number and added some clarification to the README.txt
• moved the src/recon_test/test_modelling.sh and associated input files to the recon_test_pack. It is also modified to be independent of ECAT7 and now runs again.

Other changes to tests

• expanded test_proj_data_in_memory to also test ProjDataInterfile so renamed the test to test_proj_data.

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

Major bugs fixed

• see above

Backward incompatibities

• Classes that use InterfilePDFSHeader now contain a shared_ptr<ProjDataInfo> instead of a raw pointer, removing a memory leak.
• Changes improving safety of use of shared_ptr:
  In the previous version of STIR, the use of shared_ptr allowed unsafe access to the objects (although this never happened in distributed STIR code). We now prevent this with changes to the class interface (although there is still work to do):
  • Where possible, classes that internally contained a shared_ptr<ProjDataInfo> now contain a shared_ptr<const ProjDataInfo>, and similar for DiscretisedDensity
  • get_proj_data_info_sptr used to return shared_ptr<ProjDataInfo>, but now returns shared_ptr<const ProjDataInfo>, similar for ExamInfo.
  • Corresponding constructors and some functions, including set_up, that accept shared_ptr now take a shared_ptr to a const object.
• ProjData*::copy_to and fill_from now return the updated iterator (as opposed to the size). This is like std::copy, and more convenient for reusing it.
• ModelMatrix::set_if_uncalibrated and ModelMatrix::get_if_uncalibrated renamed to ModelMatrix::set_is_uncalibrated and ModelMatrix::get_is_uncalibrated
• PlasmaData::set_if_decay_corrected and PlasmaData::get_if_decay_corrected renamed to PlasmaData::set_is_decay_corrected and PlasmaData::get_is_decay_corrected
• Remove overloaded PatlakPlot::get_model_matrix, now it can only be called using class members
• BinNormalisation::set_up() now need exam_info_sptr as input.
• BinNormalisation::apply() and undo now no longer accept start_time and end_time as they are taken from exam_info. To allow this to happen, class Bin has now an extra member time_frame(). Note that it defaults to 1. If this is incorrect, it has to be initialised explicitly (not via a constructor).
• Poisson log-likelihood hierarchy has changed for gradient methods to now use actual_compute_subset_gradient_without_penalty in the derived classes (which has an extra argument add_sensitivity). Therefore developers that created their own derived class will need to adjust their class accordingly. The public class interface is identical though.
• DetectorCoordinateMapFromfile (used by SAFIR listmode files) has been renamed to DetectorCoordinateMap and moved from listmode to buildblock.
• Scanner now has a set_up() function. This needs to be called after using any of the set_*() functions (except set_params). For the blocks/generic scanners, this will fill in the detector maps etc. You will get a run-time exception if you forgot to do this.

New functionality

• New templated functions stir::copy_to and stir::fill_from in stir/copy_fill.h which can be used to fill most STIR objects from an iterator range (or copy to). The functions normally use stir_object.begin_all() but resort to stir_object::fill_from or copy_to for a few cases where no iterators exist. We use some specialisations to try and find the fastest version.

• Introduced ProjData::standard_segment_sequence function returning 0,+1,-1,..., as used by copy_to and fill_from

• stir::read_from_file can now be called with a "node", e.g.

      auto uptr = read_from_file<VoxelsOnCartesianGrid<float> >(filename);
  

; To do this, we now need a typedef hierarchy_base_type at the top level, as otherwise we would need multiple registries. However, this was already required by write_to_file anyway.

• GeneralisedPrior and QuadraticPrior have a new method accumulate_Hessian_times_input(output, current_estimate, input). This computes the hessian of the prior at the current_image_estimate and multiplies this by input. Before this, the add_multiplication_with_approximate_Hessian(output, input) method was used that assumed a quadratic function and therefore hessian = 1.
• Exam_info, TimeFrameDefinition and PatientPosition have now a == operator
• Overloaded compute_total_ROI_values and compute_ROI_values_per_plane functions to allow the passing of aVoxelsOnCartesianGrid object instead of a Shape3D as the ROI shape. This allows the ROI volume to be constructed during preprocessing for multiple image analysis.
• KeyParser::parameter_info() now outputs vectorised keys as well
• DetectionPosition now has all comparison operators.
• LOR classes now no longer require phi at input to be between 0 and pi, nor psi to be between 0 and 2 pi. They will standardise the input automatically.

Other code changes

• store data in ProjDataInMemory in the same order as what is used by copy_to and fill_from. This enabled using straight-forward copy. (This change should not affect anyone, except if you relied on a specific order in the buffer before.)