Recent advances in iterative reconstruction for clinical SPECT/PET and CT

Abstract

Statistical iterative reconstruction is now widely used in clinical practice and has contributed to significant improvement in image quality in recent years. Although primarily used for reconstruction in emission tomography (both single photon emission computed tomography (SPECT) and positron emission tomography (PET)) there is increasing interest in also applying similar algorithms to x-ray computed tomography (CT). There is increasing complexity in the factors that are included in the reconstruction, a demonstration of the versatility of the approach. Research continues with exploration of methods for further improving reconstruction quality with effective correction for various sources of artefact.

For emission tomography (both single photon emission computed tomography (SPECT) and positron emission tomography (PET)) iterative reconstruction has virtually replaced analytical reconstruction in routine clinical practice and there is increasing interest in applying iterative techniques for reconstruction in x-ray computed tomography (CT). There is continuing research on optimising iterative reconstruction methods with particular interest in extending the current approaches to incorporate more exact, though increasingly complex models of the underlying systems. In this paper an overview of iterative reconstruction is provided with emphasis on techniques developed for application to PET and SPECT, but also with some comments on potential areas of application to CT. It is assumed that the reader is familiar with these modalities and the underlying physics. The paper focuses on the general framework for statistical iterative reconstruction rather than exploring the complete range of alternative algorithms. For additional insight the reader is referred to other published reviews on iterative reconstruction applied in nuclear medicine [1–6] or texts that provide a more extensive coverage of general reconstruction [e.g. 7].

Statistical iterative reconstruction

The most widely used reconstruction algorithm in emission tomography is the maximum likelihood expectation maximisation algorithm (ML-EM) [8,9] and its accelerated form based on ordered subsets (OS-EM) [10]. A particular advantage of these algorithms is that they are based on an underlying assumption that the measured data are subject to Poisson noise, which is true of the raw uncorrected data acquired in PET and SPECT. Alternative algorithms are more appropriate in situations where the Poisson noise assumption does not hold (e.g. see later discussion on algorithms for transmission tomography). The ML-EM algorithm involves two steps repeated iteratively: a) estimated projections are produced by forward-projection from the current object estimate for comparison with the measured projections (y_i); b) the ratio of measured to estimated projections is then back-projected and used as a multiplicative correction of the current object estimate (λ_j). This process is repeated until the estimated and measured projections are matched, resulting in a reconstruction of the activity distribution. This process is described by Equation 1:

Central to the algorithm is the system matrix (a_ij) that models the underlying physics of the detection process by describing the probabilities that photons emitted from specific object voxels will be detected at specific detector pixels. The system model can account for the detector geometry (e.g. presence of collimation) but also can easily incorporate physical factors such as the attenuation of photons and the influence of resolution or corrections such as the normalisation that accounts for variation in detector sensitivity (Figure 1). Much of the recent development work to be integrated into clinical systems has involved the extension of the system model to incorporate additional factors in order to more accurately model the emission/detection process.

Figure 1. Illustration of elements involved in iterative reconstruction. Estimated projections are established using forward projection with knowledge of the system model that can incorporate multiple parameters. The iterative update also utilises the system model but is penalised in order to favour a smooth solution; alternatively post-reconstruction smoothing can be applied.

Traditionally analytical reconstruction has been applied in transmission tomography, mainly filtered back projection (see [11,12] for description of developments in analytical and iterative reconstruction). For transmission scans the reconstruction is performed by first applying a log transform to the acquired data; hence conventional ML-EM is inappropriate as the transformed projection data are no longer Poisson distributed. Alternative iterative algorithms are therefore necessary; e.g. general methods based on use of surrogate functions [13], with specific examples being various convex algorithms [14–16]. The attraction of these algorithms is that they have properties very similar to ML-EM with improved noise characteristics and similar convergence speeds; also they are amenable to use of ordered subsets.

Modelling physical parameters

As indicated above, the system matrix should reflect as closely as possible the underlying physics of the emission and detection processes [17]. Provided the map of attenuation coefficients is known (e.g. as derived from an acquired CT scan) these values can be used to calculate the probability of detection for activity in attenuating media. The modification of the system matrix is obtained by simply weighting the probabilities by taking account of the linear attenuation from each point of emission to the detector (for SPECT) or the total attenuating path length for each line-of-response between two coincident detectors (for PET). Figure 2 illustrates, in simplified form, part of the system matrix with and without inclusion of attenuation, as well as with addition of resolution modelling. For both PET and SPECT resolution can be estimated with knowledge of the system design or directly measured; for example in the case of SPECT geometric resolution degrades linearly with distance from the collimator, in contrast PET resolution is influenced by positron range, non-colinearity of annihilation photons and detector dimensions, with, in addition, parallax effects towards the edge of the field-of-view. As more complexity is added the system matrix becomes less sparse and computational complexity, hence reconstruction time, is increased. Also convergence speed decreases with a less sparse system matrix requiring additional iterations to reach a desirable degree of convergence. Resolution modelling is now a standard option for both PET and SPECT clinical reconstruction [18–21], the main advantage being the better control of noise (due to the model being more accurate) rather than improvement in resolution, which is only achievable with sufficiently large number of iterations. In practice the noise reduction is often sacrificed in return for reduced imaging time (with less acquired counts). Note that the reconstruction algorithm (1) remains essentially unchanged for these additions to the system model.

Figure 2. Part of SPECT system matrix for cases of a simple projection model (left), model that includes attenuation (centre) and model that includes attenuation and distance-dependent resolution (right). Numbers represent the probability that emission from the source is detected in specific detector bins (projection pixels) (illustrative only).

Further modifications can be made to incorporate correction for Compton scatter in the algorithms. This requires a measurement or scatter estimation model to be incorporated in the reconstruction. A range of models have been developed for both PET and SPECT including analytical models or simplified models based on assumptions regarding the probability of scatter (see [22] for a recent review). Modelling scatter can be difficult as it depends on the distribution of both the activity and attenuation and so is object dependent. In the case of SPECT direct Monte Carlo estimation, optimised for speed of computation using convolution forced detection, has been demonstrated to be feasible in acceptable time [23]. More commonly a measured scatter estimate is incorporated, utilising narrow energy windows adjacent to the photopeak (commonly referred to as triple energy window (TEW) acquisition [24]). Whether measured or modelled, the scatter estimate can be introduced by addition of estimated scatter in the forward projection step of the ML-EM algorithm [25], either assuming a constant scatter estimate (ŝ) or refining the term in subsequent iterations (ŝ_k):

Similarly the system matrix can be easily modified to incorporate additional information such as motion parameters derived from external detection devices. The ease of incorporating the motion depends on the nature of motion and the timing of acquisition, which is quite different for PET, which uses a stationary ring of detectors, versus SPECT, which typically involves detector motion. For example in brain PET, rigid motion of the head can be directly measured at a relatively high sampling frequency and the position of the reconstructed volume can be adjusted prior to incorporation of additional counts in the amended pose (achievable by using a list-mode ML-EM algorithm [26] that accommodates individual photon updates). Clearly computational complexity increases but full list-mode reconstruction that corrects for rigid motion is certainly feasible [27,28].

Controlling noise

A further consideration in iterative reconstruction is the control of noise which increases as iteration number increases. Stopping at an early iteration is undesirable as convergence for all points in the reconstructed volume is not guaranteed. Post-reconstruction smoothing can be applied with consequential degrading of resolution although adaptive smoothing can be applied in an attempt to preserve edge resolution. The appeal in post-smoothing is the widespread availability of software based on well-understood principles, as well as the results which demonstrate a reasonable degree of position-independence for resolution.

Alternatively Bayesian algorithms can be used to introduce additional information regarding the expected noise via a ‘prior’ term. The most commonly used algorithm is the maximum a-posteriori (MAP) algorithm [29] where a ‘penalty’ term is determined, utilising the previous object estimate in a one-step-late process [30]. The cost function to be maximised is modified to include the penalty term U(λ), which is usually a function of the values in the immediate neighbourhood of each voxel. The equations used appear quite similar to ML-EM with an additional term providing a constraint to the solution that penalises differences in neighbouring values, hence reducing noise. The one-step-late MAP algorithm is described by Equation 3:

Optimising the choice of the penalty term as well as the choice of the constant (β) is a topic of continuing research (e.g. see [31]). As in post-reconstruction smoothing an adaptive term can be adopted where edges, as determined either from a separately acquired anatomical image or directly from the local distribution of counts, are preserved (see section below).

Taking advantage of new design features

The flexibility in accommodating additional information is one of the appeals of iterative reconstruction. For example in time-of-flight PET [32,33] where the difference in arrival time of the two annihilation photons is estimated, this additional information is used to amend the system matrix, since the origin of emission is more precisely known. The reconstruction is intimately linked to the development of time-of-flight technology. Likewise there is an increasing diversity in the range of instruments being designed for organ-specific imaging (e.g. compact systems for specific use in cardiac [34–36] or breast [37] imaging). Some of these systems have quite unconventional geometry and acquisition timing but the reconstruction algorithms can easily be adapted to incorporate these changes.

Similarly in development of CT for specific applications there can be severe constraints in the system design that limit the acquisition. Examples are systems designed for interventional work where complete rotation is not possible or systems designed for low-dose mammography. There is increasing awareness that sufficient conditions for reconstruction can be achieved without necessarily acquiring projections for the complete object (see [11] in the context of analytical reconstruction). Iterative reconstruction provides flexibility in accommodating unusual acquisition geometry, which encourages a wide range of possible design options (e.g. systems for tomosynthesis [38]).

Solving for additional parameters

There are instances when the goal is not only to reconstruct an unknown activity distribution but to derive further parameters or at least to allow for the existence of further variables. Motion correction based on external measurement was mentioned briefly above; however one possibility is to reconstruct with assumption of an underlying motion model, solving simultaneously for the motion parameters and the reconstructed activity [39]. Additional research aims to solve for registration, or segmentation parameters in conjunction with reconstruction (e.g. see recent work in [40,41]).

A similar approach is used to determine tracer kinetic parameters in conjunction with a time sequence of acquired images. Rather than reconstruct individual frames (noisy due to low counts) with subsequent tracer kinetic modelling, the tracer kinetic parameters can be directly determined as part of the multi-frame reconstruction process. In order to control noise a constrained solution to the time-activity curve fitting for individual voxels is achieved by assuming a limited number of basis functions can describe the kinetics rather than freely ranging parameters [42]. These approaches, though computationally demanding, offer some potential advantage in terms of signal to noise properties.

An interesting and challenging case where additional parameters are sought is the simultaneous reconstruction of an emission distribution and attenuation map from emission-only data. Since observed projections are influenced by different degrees of attenuation at different projection angles, consistency conditions can be applied when solving for the combined emission and attenuation data [43]. A promising alternative is to utilise measured scatter information to assist with the solution, based on a scatter model that is attenuation-dependent [44,45]. A computer-efficient solution would be very attractive in situations where alternative forms of attenuation measurement are not available (or in some cases not even feasible).

Incorporating information from other modalities

As briefly outlined above, noise control is an important consideration especially in reconstruction of emission tomography data. Indiscriminate smoothing results in loss of resolution and so attempts are made to constrain smoothing to sub-volumes where it can be assumed that activity is constant, based either directly on the reconstructed values or more commonly within anatomical boundaries as defined by alternative imaging modalities. Increasingly dual modality systems are available for clinical use where anatomy is directly measured by either CT or magnetic resonance imaging (MRI). The anatomical images are also used for attenuation correction and for anatomical localisation of the observed activity distribution. An example of use of anatomical information is the algorithm originally proposed by Bowsher et al. [46], where a limited number of neighbouring voxels are included in the penalty term, based on similarity of voxel values in the second modality, rather than simply including all neighbouring voxels weighted according to their distance to the central voxel. Additional weighting can be used to further support the possibility of activity variation, rather than noise, accounting for local variation in voxel values [e.g. 47] (Figure 3). Recently researchers have also explored alternative methods for incorporating information from multiple sources. For example joint entropy can be used as a cost function, an approach commonly applied for image registration, but adapted for use in image reconstruction [48–51].

Figure 3. Brain simulation illustrating (top to bottom) raw emission data; co-registered MR; conventional MLEM reconstruction; reconstructions using a conventional Bowsher prior (BP) and a modified prior with non-local weighting (NLBP) (image courtesy Daniil Kazantsev, UCL).

Potential benefits for CT

Traditionally CT has been reconstructed using filtered back projection (FBP), with both resolution and noise superior compared to emission tomography. Nevertheless there is increasing evidence that iterative reconstruction can also be useful for CT. There is potential to directly correct for relevant physical parameters such as beam hardening and scatter, including use of polyenergetic x-rays, with potential to reduce artefacts [52,53]. There is increasing interest also in improving noise structure (Figure 4), with the potential to use lower radiation exposure without degrading noise (similar to the reduced scan-time options for PET and SPECT) [54,55]. Also novel acquisition geometries are being explored for specific applications where either angular sampling is limited or truncation unavoidable; examples are in interventional systems, organ-specific designs where complete rotation is impractical or in direct dose verification during radiotherapy using cone-beam CT [e.g. 56]. In these cases iterative reconstruction tends to be more immune to artefact than analytical reconstruction. As in the case of emission tomography prior information can be incorporated (e.g. as a noise constraint) or additional information derived from additional measurement (e.g. motion [57–59]). Clearly due to the larger matrix sizes typically used in CT this places even greater demands on computational speed.

Figure 4. Coronal CT slice acquired at low dose, reconstructed with filtered back projection (left) and iterative statistical reconstruction (right) (image courtesy GE Healthcare).

Discussion

Although iterative algorithms offer gains in signal to noise and reduced artefacts there is need to recognise several limitations concerning these techniques. The convergence rate can vary with position in the reconstructed volume and is to some extent dependent on the complexity of the immediate surrounding sub-volume. Provided a sufficiently large number of iterations is used this is less problematic, however this can result in a slow and noisy reconstruction. Furthermore there can be a positive bias in areas of low activity, especially if there is correction for scatter, due to the non-negativity constraint of ML-EM reconstruction which suppresses any negative values. A further source of possible bias relates to the highly non-linear and position-dependent results that can arise when incorporating information from other modalities. For example, estimating and correcting for residual partial volume effects can prove quite difficult. As a result care needs to be taken when computing quantitative values.

As the complexity of the system model increases, computation time increases per iteration and an additional iterations are required. Consequently, until recently, reconstruction time has been the main barrier to implementing complex algorithms in routine practice. Faster computers now make these algorithms feasible, and in addition there is increasing use of fast graphics cards that provide significantly faster reconstruction than conventional processors [e.g. 60].

Meaningful evaluation and comparison of algorithms continues to be problematic as optimisation usually depends on the specific application, with subtle iteration-dependent differences in contrast and noise. Determining optimal procedures for evaluation continues to be an active area of research but detailed discussion is beyond the scope of this article (see [61]). Choosing the optimal algorithm from the diversity of available algorithms poses a significant continuing challenge for practitioners. In addition as the range of available algorithms and options increases it becomes increasingly difficult to compare and standardise results across institutions, which raises particular difficulties when designing multi-centre clinical trials.

Summary

Iterative reconstruction is now a standard method of reconstruction in clinical use, resulting in improved image quality for both emission and transmission tomography. Research continues on further extension of algorithms to account more exactly for the underlying physical characteristics of the detection systems and to account for remaining factors that produce artefacts. Research also continues on comparison and optimisation of algorithms for specific clinical applications. There is also increasing use of iterative algorithms in other areas of medical imaging (e.g. dealing with under-sampled MRI data via compressed sensing [see 62], use of prior information in optical tomography [63]); coverage, however, is beyond the scope of this article.

Acknowledgements

I acknowledge helpful input to this paper from Johan Nuyts (KUL, Louvain), Alexander Bousse (UCL) and Daniil Kazantsev (UCL). This paper was supported by EPSRC research grant number EP/G026483/1. UCL/UCLH receives a portion of research funding from the UK Department of Health's NIHR Biomedical Research Centre funding scheme.

Declaration of interest: There are no conflicts of interest to be declared.