Impact of motion compensation and partial volume correction for ¹⁸F-NaF PET/CT imaging of coronary plaque

To evaluate the performance of the MoCo and PVC methods in a scenario in which we have a full knowledge of the motion vector fields (MVFs) we used realistic simulations of PET obtained from segmentations of real MR images (Tsoumpas et al 2011a). First, an ultra-short time-echo (UTE) sequence was used in order to separate the anatomical structures of the patient and to create the attenuation maps. Second, a dynamic MR sequence (35 dynamic frames, 0.7 s duration each) acquired during normal breathing was used to estimate the simulated MVFs (Buerger et al 2012). The segmented regions of the MR images were assigned with NaF-PET standardized uptake values (SUV) as typically measured in a clinical PET/CT acquisition (Win and Aparici 2014, Oldan et al 2015), yielding a 3D emission map (table 1). Four MVFs, corresponding to different respiratory phases and obtained from the dynamic MR sequence previously mentioned, were then applied to the segmented PET image to generate a 4D segmented volume. This volume is later projected by using the STIR projector based on Siddon's ray tracing algorithm, with 20 rays per tube of response, in order to obtain the simulated gated sinograms. The number of rays per tube of response was selected to approximate resolution blurring effects. Attenuation coefficient projection data were estimated from the segmented attenuation images. The effect of scatter was simulated analytically using the single scatter simulation formula (Watson 2000) as implemented in STIR (Tsoumpas et al 2004). Random events were approximated as a uniform background, summing in total a 30% of the simulated counts. Poisson noise was included to simulate two different levels of acquired counts: 500 (high) million and 100 (low) million counts. With an injected dose of 400 MBq, 500 million counts would correspond to a 30 min scan, while 100 million counts would correspond to a 6 min scan. Further details of the process can be found in Tsoumpas et al (2011a).

For the evaluations of the two different PET systems simulated in this study, variable matrix and voxel sizes were used. Biograph TPTV PET/CT: matrix size 336  ×  336  ×  109, voxel size 2.005  ×  2.005  ×  2.031 mm³; Biograph mCT: matrix size 400  ×  400  ×  109, voxel size 2.0364  ×  2.0364  ×  2.027 mm³. STIR projector and back projector were adapted to the geometry model and 3D sinogram data storage corresponding to each system. No ToF information was used in the mCT reconstructions.

The OSEM and MoCo reconstructions were performed by using 5 iterations and 21 subsets, storing the reconstructed images after each iteration, to evaluate the convergence of the iterative algorithm with the simulated datasets. The relatively large number of iterations was chosen to increase the contrast of the small coronary plaques (Robson et al 2017). A 2 mm Full Width Half Maximum (FWHM) Gaussian filtering was applied after each iteration in all reconstructions, following as reference the settings used in other similar studies (Rubeaux et al 2016, Robson et al 2017). No regularization was included in any of the static or MoCo reconstructions.

A spherical shape of the plaque, with 2 mm radius (33.5 mm³ volume), was simulated in all cases. The assigned SUV values for the calcified plaque were set to obtain lesion-to-background ratios (LBR) of 10, 20, 50 and 70 (background SUV  =  0.8), following the values suggested in Delso et al (2011). In order to evaluate the variability of the results, five simulations with different seed numbers were analyzed. The seed numbers were used to create different noise realisations when adding Poisson noise to the simulated gated-sinograms.

The figures of merit used for the evaluation of the simulated reconstructions were: LBR using the maximum value within the segmented lesion (LBR_max), LBR using the mean activity within the segmented lesion (LBR_mean), activity bias (%) in each lesion and image noise (%), measured as the standard deviation in a uniform region. The activity bias within the lesions was computed as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm bias}\left( \% \right)=~\frac{\mathop{\sum }_{v\in {{L}_{i}}}\frac{\left( {\rm measure}{{{\rm d}}_{v}}-~{\rm re}{{{\rm f}}_{v}} \right)}{{\rm re}{{{\rm f}}_{v}}}}{N{\rm voxel}{{{\rm s}}_{Li}}}\cdot 100\nonumber \end{align} \tag{ 1 }$$

where measured_v is the activity measured in each voxel within the segmented lesion L_i, ref_v is the activity in the reference image, down-sampled to the PET resolution, and Nvoxels_Li is the total number of voxels within the segmented lesion.

The errors for LBR_max were calculated as the standard deviation of the LBR_max values obtained in five independent case-based simulations, weighted with the noise measured in the uniform region. The errors on the LBR_mean, bias and noise values were computed as the standard deviation obtained from the five independent simulations.

To evaluate the effect of the MoCo  +  PVC methods on the quantification of the PET images, we calculated the relative changes in the LBRs (ΔLBR) after applying MoCo and MoCo  +  PVC. The ΔLBR was calculated as:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta {\rm LBR}\left( \% \right)=\frac{{\rm LB}{{{\rm R}}_{{\rm corrected}}}-{\rm LB}{{{\rm R}}_{{\rm static}}}}{{\rm LB}{{{\rm R}}_{{\rm static}}}}\cdot 100\nonumber \end{align} \tag{ 2 }$$

where LBR_corrected is the LBR after applying the MoCo or MoCo  +  PVC methods and LBR_static is the LBR measured in the static image reconstructed without any corrections.

2.2.1. Phantom experiments

Acquisitions of a human thorax phantom (Fieseler et al 2013), capable of mimicking respiratory and cardiac motion, were also used for the evaluation of the MoCo  +  PVC approach. This thorax phantom contained three spherical plaque-type lesions of 36, 31 and 18 mm³ inserted (figure 1). The background activity was set to 4 kBq ml⁻¹, while the plaque-type lesions were filled with a LBR of 70:1. No increased myocardial uptake was simulated in any of the experiments, in order to mimic the conditions of clinical imaging with the NaF radiotracer. The respiratory motion range was set to 2 cm, while the cardiac motion was set to a maximum change of 60 ml blood volume. The following acquisitions were performed:

• Exp1. Static acquisition with no motion, 10 min listmode. Systole and expiration. 112 million counts.
• Exp2. Only cardiac motion, 10 min listmode. Cardiac cycle 1.2 s (0.6 s systole, 0.6 s diastole). 124.1 million counts.
• Exp3. Respiratory and cardiac motion, 16 min listmode. Respiratory cycle 4.6 s (1.6 s inspiration, 3.0 s expiration). Cardiac cycle 1.2 s (0.6 s systole, 0.6 s diastole). 163.5 million counts.

Zoom In Zoom Out Reset image size

All experiments were conducted using a Biograph mCT PET/CT system, located at Münster University (Germany). Images were reconstructed following the same settings as in the simulated datasets. In Exp1, static OSEM and static OSEM  +  PVC reconstructions were performed. This experiment is used as reference standard for the accuracy measurement of the motion estimation and compensation techniques implemented in this work. Cardiac motion was modelled in Exp2, while both respiratory and cardiac motions were modelled in Exp3 (see section 2.4 for further details). LBR_max and LBR_mean values were evaluated for all the lesions.

2.2.2. Patient studies

A single respiratory (R-MoCo) or cardiac (C-MoCo) was used in the simulations and in phantom experiment 'Exp2', respectively. We used the STIR C++ reconstruction software (Thielemans et al 2012). STIR is open source software for use in tomographic imaging, which is focused especially on iterative image reconstruction in PET and SPECT. STIR supports off-line image reconstruction by using several different algorithms and system geometries, it also incorporates most of the corrections traditionally performed in PET (attenuation, scatter, normalization, dead time, decay, etc), and for this reason has received widespread acceptance in the field of nuclear medicine.

Figure 2 shows a general scheme of the MoCo  +  PVC reconstruction framework implemented in this work. For the MoCo reconstruction, we use the emission and attenuation sinograms (${{Y}_{bg}}$ and A_bg) and the MVFs (${{\widehat{W}}_{g}}$ and $\widehat{W}_{g}^{-1}$) for each gate g. The MVFs are used to generate a series of gated attenuation images from the attenuation image that corresponds to the reference gate. These images are then projected into sinogram space, thus, obtaining the gated attenuation sinograms used in MCIR. Finally, the PVC is performed for the MoCo images, obtaining a MoCo  +  PVC image (see section 2.5).

Zoom In Zoom Out Reset image size

For the motion compensation we make use of an iterative algorithm with a common regularization approach (equation (1)): the ordered subsets maximum a posteriori one-step-late, i.e. OSMAPOSL (Green 1990a, 1990b), which is freely available within the STIR library (Thielemans et al 2012). This algorithm was extended (Tsoumpas et al 2013) to include motion information in order to reconstruct all respiratory or cardiac gates into one reference frame. The corresponding MoCo-IR equation for incorporating the motion compensation within the system matrix is (Tsoumpas et al 2013):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \begin{array}{@{}lllllllllllllllllllllll@{}} & \Lambda _{v}^{\left(s+1 \right)}=~\Lambda _{v}^{(s)}\frac{1}{\mathop{\sum }_{b\in {{S}_{l}},g}\mathop{\sum }_{{{v}^{\prime}}}\widehat{W}_{{{v}^{\prime}}g\to v}^{-1}{{P}_{{{v}^{\prime}}b}}{{A}_{bg~}}}\nonumber \\ &\qquad\qquad\,\,\,\times \underset{b\in {{S}_{l,g}}}{\mathop \sum }\,\underset{{{v}^{\prime}}}{\mathop \sum }\,\left(\widehat{W}_{{{v}^{\prime}}g\to v}^{-1}~{{P}_{{{v}^{\prime}}b}}~\frac{{{Y}_{bg}}}{\mathop{\sum }_{\widetilde{v}}{{P}_{b\widetilde{v}}}~\mathop{\sum }_{\widetilde{v}}{{\widehat{W}}_{\widetilde{v}\to \widetilde{v}g}}\Lambda _{\widetilde{v}}^{(s)}~+~\frac{{{B}_{bg}}}{{{A}_{bg}}}} \right) \nonumber \\ \end{array}\nonumber \end{align} \tag{ 3 }$$

where $\Lambda _{v}^{(s)}$ is the estimated radioactivity at voxel v and sub-iteration s. ${{Y}_{bg}}$ is the number of measured coincident photons of each detector pair b that belongs to the lth subset S and gate g. S_l is the lth subset of the projection space, which is divided into L total subsets. s is the sub-iteration number. P_bv is the system projection matrix. $\widehat{W}$ and ${{\widehat{W}}^{-1}}$ represent the forward and backward warping operations of the image that move the activity form one location to another. The operator ${{\widehat{W}}^{-1}}$ acts to warp all respiratory or cardiac gates to the reference gate, while the operator $\widehat{W}$ 'unwarps' the reference gate to each gate. Finally, A_bg and B_bg are the attenuation coefficient and background terms for each bin and gate. No regularization was used in any of the reconstructions in order to compare images reconstructed using exactly the same settings. The above mentioned method for motion compensation is expected to work well, provided that accurate models for the respiratory and/or cardiac motion are available.

The MVFs used for the single respiratory (R-) or cardiac (C-) MoCo reconstruction were estimated from the co-registration of previously reconstructed gated PET images to a reference gate (end-expiration for respiratory gating and diastole for cardiac gating). The MIRT (Medical Image Registration Toolbox) software (Myronenko and Song 2010) was used in Matlab for the co-registration of the gated images.

In the simulated datasets (see section 2.1), the accuracy of the estimated MVFs at different levels of noise (MVF obtained from ideal gated images with no noise and gated images obtained from acquisitions with 500 M and 100 M counts) was evaluated by direct comparison, at different volumes of interest (VOIs), of these MVFs with the original MVFs used for creating the reference gated images. Six VOIs were used for the evaluations, as depicted in figure 3.

Zoom In Zoom Out Reset image size

For the experimental data in which both respiratory and cardiac motions were present (phantom experiment 'Exp3' and all patient scans), respiratory and cardiac motion compensation (RC-MoCo) is required. In these cases we used a joint RC-MoCo, which first models the respiratory motion using markerless motion detection, and corrects for it in the projection space prior to image reconstruction, and then applies the cardiac MoCo during the reconstruction process. The proposed method is illustrated in figure 4. Given the main potential application of the method being cardiac PET, we assume that cardiac triggers are available from an electrocardiogram (ECG). In contrast, no respiratory triggers where acquired, which is in fact the case in all cardiac studies evaluated in this work. The procedure is described in the following:

• The respiratory signal is obtained by using a data-driven motion detection (Büther et al 2009), which measures the changes in count rates within a volume of interest (VOI) clearly affected by motion (the diaphragm was used in this work). We used a center-of-mass (CoM) analysis of the count-rate of the sinogram bins connected with the projection of the segmented diaphragm. The count-rate analysis is performed using a 100 ms time-division of the LM file (Lassen et al 2015).
• The respiratory signal is used for a subsequent amplitude-based (5-95-percentile of the respiratory range in the patient data, full respiratory range for the phantom studies) binning of the LM data. The respiratory gated sinograms (4 in this study) are reconstructed using the STIR software, without attenuation correction.
• The ideal projections (without noise) of the gated non-attenuation-corrected (NAC) PET images obtained in the previous step are used for the definition of motion vector fields (sMVFs) in the projection space for each gated sinogram. This co-registration was performed using the MIRT toolbox (Myronenko and Song 2010).
• These MVF are used for MoCo of each event in the listmode data. The transformed data is then stored in a new 'respiratory motion-free' listmode (LM) file. To ensure correct normalization of the new LM data, we applied a correction factor to each event in the LM file, computed as the ratio of the normalization factors of the original and the motion-compensated positions.
• The original cardiac triggers are also stored in the new LM-data and used to create cardiac-gated sinograms.
• These sinograms are then reconstructed using the MoCo method described in section 2.3. The cardiac iMVFs were obtained from a non-rigid co-registration of prior reconstructed cardiac-gated PET images, by using the MIRT software (Myronenko and Song 2010).

Zoom In Zoom Out Reset image size

Further details about the markerless respiratory motion detection and compensation can be found in Lassen (2017).

Here, we employ the LP method for PVC of the PET images (Moore et al 2012, Cal-González et al 2015). The LP method is used as a post-processing step after reconstruction. This method employs a segmentation of typically a few tissues within a small volume of interest centred on the lesion to be analyzed, as shown in figure 5 (for J  =  2 tissues). The measured emission projection counts, λ, are modelled as the sum of the projection counts from each segmented tissue (with unit-activity concentration), scaled by their respective activity, plus the counts coming from the global background outside the VOI, as shown in equation (4):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\lambda }_{i}}\equiv \underset{j=1}{\overset{J}{\mathop \sum }}\,{{A}_{j}}{{P}_{ij}}+~{{g}_{{\rm out},i}}\nonumber \end{align} \tag{ 4 }$$

where λ_i are the expected counts per sinogram bin, i, A_j the activity for each segmented tissue j within the VOI, P_ij is the resolution-blurred tissue shape function for tissue j and sinogram bin i and ${{g}_{{\rm out},i}}$ represents the background counts coming from the region outside the VOI.

Zoom In Zoom Out Reset image size

The joint likelihood of measuring a given projection data set, n, is given by the product of the Poisson probability density function for each measured projection ray, i, over all rays traversing the VOI:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle L~=\underset{i~\in ~{\rm VOI}}{\mathop \prod }\,\left[ \frac{{{{\rm e}}^{-{{\lambda }_{i}}}}\lambda _{i}^{{{n}_{i}}}}{{{n}_{i}}!} \right].\nonumber \end{align} \tag{ 5 }$$

Where ${{n}_{i}}$ are the measured counts per LOR. Taking the logarithm of the likelihood, we obtain:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm ln}\left(L \right)=~\underset{i}{\mathop \sum }\,{{n}_{i}}{\rm ln}\left({{\lambda }_{i}} \right)-~{{\lambda }_{i}}-{\rm ln}\left({{n}_{i}}! \right).\nonumber \end{align} \tag{ 6 }$$

When substituting the λ_i values with equation (4):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm ln}\left(L \right)=~\underset{i}{\mathop \sum }\,{{n}_{i}}{\rm ln}\left(\underset{j=1}{\overset{J}{\mathop \sum }}\,{{A}_{j}}{{P}_{ij}}+~{{g}_{{\rm out},i}} \right)-~\left(\underset{j=1}{\overset{J}{\mathop \sum }}\,{{A}_{j}}{{P}_{ij}}+~{{g}_{{\rm out},i}} \right)-{\rm ln}\left({{n}_{i}}! \right).\nonumber \end{align} \tag{ 7 }$$

The tissue-activities A_j are determined by maximizing the log likelihood (equation (4)) for the expected value λ_i. This can be done by taking derivatives of ln(L) with respect of each tissue-activity A_j and setting the results to zero:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle A_{j}^{(k)}=\underset{{{j}^{\prime}}=1}{\overset{J}{\mathop \sum }}\,{{\left[ \underset{i}{\mathop \sum }\,\frac{{{P}_{i{{j}^{\prime}}}}{{P}_{ij}}}{\lambda _{i}^{(k)}} \right]}^{-1}}~\underset{i}{\mathop \sum }\,\frac{{{P}_{ij}}\left({{n}_{i}}-g_{{\rm out},i}^{(k)} \right)}{\lambda _{i}^{(k)}}.\nonumber \end{align} \tag{ 8 }$$

Note that the equations in (8) cannot be solved analytically for the unknown A_j'. Instead, an iterative solution for these equations is sought, thus the values of ${{\lambda }_{i}}$, ${{g}_{{\rm out},i}}$ and A_j are updated per iteration k.

The tissue-activities obtained by the LP method are being used as a prior in a new reconstruction procedure, performed with the STIR software (Thielemans et al 2012), in order to obtain a locally PVC-PET image (figure 5). First, the system matrix values (P_ij) for each tissue and the global background g_out,i that affects the VOI are computed. Later, the LP tissue-activities are obtained and the activity in each voxel within the VOI is substituted with the tissue-activities obtained with the LP method. The resulting image is forward projected using the STIR projector and the obtained sinogram is reconstructed using STIR resulting in a PVC image (see Cal-Gonzalez et al (2017) for details).

Table 2 and figure 6 show the estimated mean MVFs values (in mm) from respiratory-gated PET images, reconstructed without noise and with the noise corresponding to a total of 500 M and 100 M counts. The reference MVFs, used for creating the simulated datasets, are also shown in the table and figure. The accuracy of the estimated MVFs, and therefore, of the R-MoCo approach, is directly related with the amount of noise in the gated PET images used for the estimation of the MVFs.

Zoom In Zoom Out Reset image size

Figure 7 shows the Static, R-MoCo and R-MoCo  +  PVC images in the Biograph mCT system. Similar images were obtained for the Biograph TPTV system. In general, we observe a significant improvement of the detectability and quantification of the small plaque lesions when the MoCo is applied. Further improvements are observed when R-MoCo  +  PVC are applied. This is demonstrated in figures 8 and 9, where we plot the LBR_max (figure 8) and bias values (figure 9) for all the plaque lesions. In addition, figure 10 shows the noise in the reconstructed images, measured from a VOI located in a region with uniform activity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 3 presents the LBR_max and LBR_mean values for the simulated plaque in each case, while their relative changes (ΔLBR) after applying the R-MoCo and PVC corrections are shown in table 4. The values in the figures and in the table demonstrate clearly the impact of motion compensation and PVC in the NaF uptake in coronary arteries, with ΔLBR_max values between 120 and 945% for the MoCo  +  PVC images and between 200 and 1119% for the tissue-activities obtained with the LP method.

The markerless detection of respiratory motion was evaluated on the moving phantom (Exp3), using the Anzai respiratory belt as reference (Anzai Medical, Japan). Figure 11(A) shows the frequency spectra obtained by means of a FFT analysis and the filtered respiratory signal, together with the signal obtained with the Anzai respiratory belt. In addition, figure 11(B) shows the frequency spectra and the filtered respiratory signal for one of the ¹⁸F-NaF scans (patient 1). As expected, we obtained a strong respiratory peak in the Fourier spectrum for the phantom experiment (regular periodic motion), while for the patient we obtained a broad peak between 0.1 and 0.4 Hz, due to the more irregular motion of the patients.

Zoom In Zoom Out Reset image size

The measured respiratory translations of the diaphragm (axial direction only) from the respiratory-gated NAC-PET reconstructions are shown in table 5, for the phantom experiment (Exp3) and for all the patient scans.

Figure 12 shows the images of the phantom with not motion (Exp1), with and without PVC; together with the static, MoCo and MoCo  +  PVC images of the moving phantom, for the acquisitions Exp2 and Exp3. The LBR_max and LBR_mean values for the three lesions are presented in table 6, while their relative changes (ΔLBR) after applying the MoCo and PVC corrections are shown in table 7.

Zoom In Zoom Out Reset image size

As expected, a significant improvement in the visual aspect of the plaques (figure 11) was achieved when MoCo  +  PVC were performed, giving images closer to the ones obtained in the static acquisition (Exp1). This assumption is reinforced with the increase observed in the LBR_max values of the plaques for the MoCo  +  PVC images. Further improvement of the LBR values can be observed when using the tissue-activities from the LP method, applied from the MoCo image (right column on table 6).

As shown in table 7, higher relative changes in the LBR values were observed when applying the RC-MoCo  +  PVC in the moving phantom with both respiratory and cardiac motions incorporated (ΔLBR_max values in Exp3 between 157 and 442% for the RC-MoCo  +  PVC images and between 381 and 614% for the LP method).

In total, 6 plaque lesions with increased ¹⁸F-NaF uptake were observed. The static OSEM, R-MoCo, RC-MoCo and RC-MoCo  +  PVC images obtained with the STIR software are shown in figure 13, together with the AC-CT image. Notably, the calcified plaques are easier to see when MoCo and PVC are applied. The PVC further improves the delineation of the lesion on PET with the lesion detected in the Ca-score CT.

Zoom In Zoom Out Reset image size

The LBR_max and LBR_mean values, obtained from the images (Static, R-MoCo, RC-MoCo and RC-MoCo  +  PVC) and from the LP method (for the RC-MoCo images) are shown in table 8. These values were computed by dividing the LP tissue-activity for the lesion and the mean activity of the background VOI. Table 9 shows the ΔLBR values when applying R-MoCo, RC-MoCo, RC-MoCo  +  PVC and the tissue-activities values directly obtained from the LP method. As expected, a significant increase in the LBR can be observed when RC-MoCo and PVC are incorporated into the reconstruction (ΔLBR_max values between 31 and 211% for the RC-MoCo  +  PVC images and between 46 and 373% for the LP method). Furthermore, we can observe a bigger effect of the PVC (higher ΔLBR values when the size of the segmented plaque is smaller).

The main limitation of the study is the lack of a proper evaluation of the RC-MoCo  +  PVC approaches with simulated data, in order to fully explore the accuracy of the estimated respiratory and cardiac MVFs and the performance of the LM-based respiratory MoCo. Such evaluation is currently under investigation and will be presented in future work.

A potential limitation of the phantom and patient studies presented in this paper is the possible lack of accuracy of the derived motion fields, due to the noisy nature of the gated PET images and intra-gate motion not taken into account. As a consequence, the increase in the LBR of the calcified lesions in the images is not as high as expected from the simulations.

Another important consideration is the number of gates chosen for the estimation of the MVFs (4 gates for both respiratory and cardiac motion modelling), which is below the standard number of gates used in clinical routine (usually 8 gates) and might have a negative impact on the performance of the MoCo  +  PVC due to uncorrected intra-frame motion. However, due to the high level of noise in our gated images, and using the results obtained from the simulations as reference, we considered that 4 gates would be the best choice to achieve a good compromise between intra-frame motion and noise. The optimal number of gates is a topic of further investigation.

In addition, the lack of time of flight (ToF) modelling in our reconstructions can be also considered a limitation of the study. This is due to the fact that using the TOF information in the reconstruction will result in better image contrast and lower noise, which can help to obtain more accurate MVFs from the noisy, gated PET reconstructions. Furthermore, the ToF information can also be used to improve a joint motion estimation and image reconstruction of PET data, as recently discussed by Bousse and colleagues (Bousse et al 2016).

Finally, correction for non-periodic motion will most likely require the generation of adaptable respiratory and cardiac motion models that will be integrated with list-mode reconstruction (King et al 2011).