Quantitative image reconstruction for total-body PET imaging using the 2-meter long EXPLORER scanner

The simulated EXPLORER system configuration and simulation parameters are listed in table 1. The scanner consists of 36 axial block rings with an axial gap of 3.42 mm (one crystal pitch) between adjacent rings. Each ring has 48 $51\times 51$ mm² detector modules forming a ring of 800 mm in diameter. Each detector module consists of an array of $15\times 15$ LSO crystals with a crystal pitch of 3.42 mm. The time-of-flight (TOF) resolution of 530 ps was chosen to mimic the timing performance in a Siemens Biograph mCT PET scanner.

Table 1. System configuration of the simulated EXPLORER scanner.

System parameters	Value
Scintillator material	LSO
Crystal size (mm)	$3.34\times 3.34\times 20$
Crystal pitch (mm) (+80 μm reflector)	3.42
Number of crystals per block detector	$15\times 15$
Number of block detectors per ring	48
Number of block rings	36
Ring diameter (mm)	800
Transaxial FOV (mm)	700
Gap size between adjacent block rings (mm)	3.42 (one crystal pitch)
Axial FOV (mm)	1966
Energy resolution	13%
Energy window (keV)	[435, 650]
TOF resolution (ps)	530
Timing bin size (ps)	25
Coincidence timing window (ns)	±5.5

Statistically-based iterative image reconstruction methods can provide improved image quality by using accurate statistical and physical models of the PET imaging process. In general PET data $\boldsymbol{y}\in {{\mathcal{R}}^{M\times 1}}$ can be modeled as a collection of independent Poisson random variables with the expectation, $\bar{\boldsymbol{y}}\in {{\mathcal{R}}^{M\times 1}}$, related to the unknown radiotracer distribution $\boldsymbol{x}\in {{\mathcal{R}}^{N\times 1}}$ through an affine transformation with a system matrix $\boldsymbol{P}\in {{\mathcal{R}}^{M\times N}}$ (Qi and Leahy 2006)

$$\begin{eqnarray}&&\bar{\boldsymbol{y}}=E\left[\,y|x\right]=Px+\bar{s}+\bar{r}\end{eqnarray} \tag{ 1 }$$

The (i, j)th element of $\boldsymbol{P}$, p_ij, models the probability of a radioactive decay from the jth image voxel being detected by the ith LOR. $\bar{\boldsymbol{s}}$ and $\bar{\boldsymbol{r}}$ denote the expectations of scatters and randoms, respectively.

Since the EXPLORER scanner has more than 50 billion LORs and records time-of-flight information, list-mode data format is preferred for efficient data processing. An additional consideration for EXPLORER is that the axial resolution can be degraded by the parallax effect between crystal pairs with a large ring difference as illustrated in figure 2. A more oblique LOR (red color) penetrates more adjacent crystals than a LOR of less ring difference (blue color). It is similar to the radial blurring effect at large radial offset in a transaxial plane. Figure 2 (lower right) shows reconstructions of a simulated point source without any resolution modeling with a maximum block-ring difference (MBRD) of  ∼30, 110 and 200 cm, respectively. Clearly we observe worse axial resolution with longer tails as we increase the maximum ring difference, although the resolution degradation effect is less dramatic because the effect of oblique LORs is mitigated by direct LORs passing through the same point. Therefore, it is necessary to model the resolution loss in both radial and axial directions in the list-mode image reconstruction. We adopt an image domain resolution model by factoring the system matrix $\boldsymbol{P}$ into three major components (Zhou and Qi 2014):

$$\begin{eqnarray}&&\boldsymbol{P}=\boldsymbol{D}\boldsymbol{G}\boldsymbol{R}\end{eqnarray} \tag{ 2 }$$

where $D=\text{diag}\{n\}$ is a diagonal matrix containing normalization and attenuation factors, $G$ is a simplified geometric projection matrix calculated by a ray-tracing algorithm (Siddon 1985) with a TOF kernel, and $R$ is an image blurring matrix which models the resolution degradation effects, such as the positron range, photon acollinearity, and detector responses including inter-crystal penetration and inter-crystal scatter effects.

Zoom In Zoom Out Reset image size

For simplicity, the image domain point spread functions (PSFs) are considered to be separable in the transaxial plane and axial direction, i.e., the PSF at each point satisfies

$$\begin{eqnarray}&&\text{PSF}(x,y,z)=\text{PS}{{\text{F}}_{t}}(x,y)\text{PS}{{\text{F}}_{a}}(z).\end{eqnarray} \tag{ 3 }$$

Note that both the transaxial PSF and axial PSF are spatially variant. However, we assumed that the $\text{PS}{{\text{F}}_{t}}$ is rotationally invariant, i.e., two points with the same radial distance to the FOV center share the same transaxial PSF after rotation, and the $\text{PS}{{\text{F}}_{a}}(z)$ is kept invariant within the same axial plane. Therefore, we only need to estimate the transaxial PSF along a radial line and the axial PSF along the scanner z-axis. We simulated 101 point sources along the radial direction and 265 point sources along the axial direction (¹⁸F source in water for modeling the positron range) with a spacing of 3.42 mm which is equal to the voxel size used in reconstructions. We used a list-mode TOF ML-EM (LM TOF ML-EM) algorithm to reconstruct the 366 point sources individually, and then estimated the separable blurring kernels in the whole FOV by fitting a 2D Gaussian mixture model (two Gaussians) in the transaxial plane and 1D Gaussian model in the axial direction, respectively. Figure 3 shows examples of the estimated blurring kernels and the FWHM (full width at half maximum) and FWTM (full width at tenth maximum) values are listed in table 2 for selected radial locations. The total storage size of the PSFs is less than 50 MB for the EXPLORER.

Zoom In Zoom Out Reset image size

Normalization is required in reconstruction to compensate for the variations in detector efficiencies and geometric factors that are not modeled in the geometric projection matrix. To obtain the normalization factors, we simulated a uniform cylinder with a dimension of 70 cm in diameter and 200 cm in length using SimSET. To reduce the simulation time, we used a $1\times 1\times 1$ voxel for the emission map and attenuation map and set the parameter of 'Target Cylinder' to the above-mentioned dimensions to constrain the source inside the cylinder. ¹⁸F source in 'air' was used without adjustment for positron range and photon collinearity. After  ∼2800 cpu-hours of simulation time, a total of 12.8 billion true coincidence events were acquired. To increase the statistics of the normalization factor in each individual LOR, we took advantage of the geometric symmetry in the scanner. For each LOR, there is a 48-fold rotational symmetry in the transaxial plane by rotating the LOR over $M\times 360/48$ degrees (M  =  1,.., 47, and 48 is the number of transaxial blocks) and another 2-fold reflection symmetry for the LORs belongs to the same block pair in the transaxial plane. In addition, there is a $N\left(=36-\text{BRD}\right)$-fold parallel symmetry in the axial direction, where BRD is the block ring difference ranging from 0 to 35. Therefore, the counting statistics of each LOR can be increased by $96\times N$-fold. Figure 4 shows an example of the LORs sharing the same normalization factor, where the red line represents the original LOR, the blue, green and purple ones represent the symmetric LORs. Figure 5(a) shows the original sinogram of the cylinder with axial symmetry only, which has no obvious pattern because of the low counts per pixel. Figure 5(b) shows the sinogram after averaging by both the axial symmetry and transaxial symmetry. The pronounced block pattern was caused by the inter-crystal scatters and the energy threshold, which effectively reduces the detection efficiency of edge crystals in a block. The final normalization factors were computed as the ratio between the averaged SimSET sinogram and the forward projection of the same cylinder using $\boldsymbol{G}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Due to the long axial FOV of EXPLORER, photons may have longer paths through scattering material, resulting in greater opportunities for multiple scattering than may be seen in conventional PET scanners. To model multiple scatters, we used a Monte Carlo based scatter estimation in this study. We assume that the scatter sinogram is relatively smooth after proper correction by detector efficiencies. Therefore, we estimate the expectation of the scatter sinogram by histogramming the scatter events in Monte Carlo simulation by detector block pairs and also using a coarse TOF bin size (e.g. 225 ps here). This effectively increases the counts by ${{15}^{4}}\times 225/25\approx 4.5\times {{10}^{5}}$ fold, resulting in a nearly noise-free scatter sinogram. For the simulated EXPLORER with 36 block rings and 48 detector modules per ring, the storage size of the 5D³ block-pair based scatter sinogram is 197 MB in floating point. In the list-mode image reconstruction, we computed the scatter mean for each event using a 5D linear interpolation between the adjacent block-pairs and TOF bins, and then scaled to the timing bin (e.g. $25~\text{ps}$ here) for reconstruction.

For random correction, we use the following formula to calculate the mean random rate (cps) (Knoll 2000)

$$\begin{eqnarray}&&{{R}_{ij}}=2\tau {{S}_{i}}{{S}_{j}}\end{eqnarray} \tag{ 4 }$$

where τ is the coincidence timing window, and S_i and S_j are the singles rates of crystals i and j, respectively. For TOF data, the formula is changed to

$$\begin{eqnarray}&&R_{ij}^{TOF}={{\tau}_{\Delta}}{{S}_{i}}{{S}_{j}}\end{eqnarray} \tag{ 5 }$$

where ${{\tau}_{\Delta}}$ is the TOF bin size (25 ps in our simulated study).

To evaluate the system performance, we conducted computer simulations using the SimSET Monte–Carlo toolkit (Harrison et al 2006). We employed the XCAT 2.0 anthropomorphic adult male phantom (Segars et al 2010) (figure 6) with a height of 178 cm in the simulation. To evaluate the capability of low dose imaging provided by EXPLORER, we modeled an injected activity of 25 MBq (∼675 μCi), which is about 1/20th of the standard dose administered in clinical PET scan protocols. To simulate a 20 min ¹⁸F-FDG PET scan starting 60 min post-injection, different time-activity curves were generated for 14 major organs and tissues (Zhang et al 2014a). In addition, two spherical lesions of 10 mm in diameter were simulated inside of the liver and lung with a local contrast of 3.2 and 14.4, respectively.

Zoom In Zoom Out Reset image size

To generate random coincidences, we ran a separate SimSET simulation in the 'SPECT' mode and then doubled the event rates to get the coincidence singles rates for PET. The ¹⁷⁶Lu source in LSO cannot be simulated directly using SimSET. Instead we used a measured singles rate of 105 cps cc⁻¹ (within the energy window of 435–650 keV) based on Siemens mCT block detectors. Therefore, the singles rate from LSO background is about 24 cps per crystal (volume of 0.23 cm³), which was added to the singles events from the object. Then we generated a Poisson random variable with the mean of ${{R}_{ij}}=2\tau {{S}_{i}}{{S}_{j}} \Delta T$ (Knoll 2000) to simulate random coincidences in each LOR, where $\Delta T$ is the scan duration. The maximum ring difference was set to 20 axial blocks (115 cm) to maximize the noise equivalent count rate (NECR) based on our previous studies (Poon et al 2012, Poon 2013). A constant coincidence timing window of $\tau =5.5$ ns was used to cover the entire FOV. In practice, a ring-difference dependent variable timing window can be used to reduce the number of random events, but the fraction of randoms within the FOV and the image quality will not be affected when TOF information is used.

Table 3 lists the three imaging protocols that were evaluated to compare the performance of the EXPLORER for low-dose static imaging with a clinical scanner. The first imaging protocol considered a 4 block-ring PET scanner at a single bed position that recorded data for 20 min over the liver region (4BR-single, figure 7(a)). The second imaging protocol was a whole-body multi-bed scan using the 4 block-ring system (4BR-multi, figure 7(b)). Each bed scan lasted 0.6 min and the total imaging time was also 20 min. This protocol is similar to whole-body imaging with continuous bed motion. The third imaging protocol was using the EXPLORER (36BR, figure 7(c)) and also corresponded to a 20 min scan. The first two sets of simulation data (4BR-single and 4BR-multi) were generated by sorting the events from the EXPLORER simulation with a maximum block-ring difference of 3 (MBRD  =  3), which mimics a current clinical scanner with an axial FOV of 22 cm. For each imaging scenario, we estimated the respective image domain resolution model using point source scans. Ten independent, identically distributed noisy realizations were generated for each protocol to estimate the uncertainty in the region of interest (ROI) quantification.

Zoom In Zoom Out Reset image size

To evaluate the image quality, we calculate the contrast recovery coefficients (CRC) of the liver and lung lesions and the background variability (STD) from the reconstructed images. The true lesion masks were used as the lesion ROIs. The background ROIs were drawn in the liver and lung, respectively, excluding the lesion and boundary voxels. The CRC was calculated by

$$\begin{eqnarray}&&\text{CRC}\left(\% \right)=\frac{\left(\text{lesion}~\text{mean}\right)/\left(\text{background}~\text{mean}\right)-1}{\text{true}~\text{contrast}-1}\times 100 \% \end{eqnarray} \tag{ 6 }$$

and the STD was defined as the standard deviation of the voxel values inside the background ROI.

The software of the 3D TOF list-mode OS-EM algorithm was coded in C++ and compiled using the GNU GCC (g++) compiler with the option 'O3' for optimization. Ten iterations and 10 subsets were used for all image reconstructions which is enough to achieve effective convergence (see results in the next section). The image dimension were $195\times 195\times 527$ with a voxel size of $3.42\times 3.42\times 3.42$ mm³. Processing 3.17 billion coincidences took 8.7 min per iteration (6.1 million events per second) on a computer with dual sixteen-core CPUs (Intel(R) Xeon(R) CPU E5-2683 v4 @ 2.1 GHz).

We first plot in figure 8(a) the numbers of detected true events and scattered events as a function of ring difference for the anthropomorphic phantom simulation. Clearly both numbers decrease rapidly with increasing ring difference due to the fact that there are fewer number of sinograms and stronger object attenuation at larger ring differences. The small peaks in the curves are caused by the axial gaps between block rings. Figure 8(b) shows a trend of increasing scatter fraction in individual sinograms as a function of the ring difference (red curve), but the overall scatter fraction in the fully 3D PET data remains relatively flat across all maximum ring differences (blue curve). The slight decrease in scatter fraction for ring difference less than 25 cm is because the anthropomorphic phantom is only 178 cm tall and does not cover the whole axial FOV of the scanner.

Zoom In Zoom Out Reset image size

Furthermore we examined the contributions of multiple scatters in the EXPLORER. Figure 9 plots the percentages of single scatters, double scatters, triple scatters and other scatters as a function of ring difference. Figure 10 shows the axially summed scatter sinograms. Here the double and triple scatters are defined by the total number of scattering interactions undergone by the two coincident photons inside the object. For example, a double scatter means either one photon is scattered twice or both photons are scattered once each. We can see that the fraction of multiple scatters increases as a function of ring difference with the percentage of single scatters decreasing from 85% in direct planes to 75% at the largest ring difference.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To verify the accuracy of the 5D linear interpolation, we compared the axially summed scatter sinogram and estimated scatter mean from the 5D interpolation in figure 11. The difference sinogram between the estimated scatter mean and MC scatters is shown in figure 11(c) and the profile comparison in figure 11(d). It shows that the block-pair based average with 5D interpolation provides a good estimate of the scatter distribution.

Zoom In Zoom Out Reset image size

Figure 12(a) shows the singles rate map from the SimSET simulation in combination with the LSO background, which represents S_i and S_j in equation (5). The horizontal and vertical axes correspond to the transaxial and axial crystal indices, respectively. Figure 12(b) shows the axial sum of the estimated mean of the random events, which is relatively uniform across the FOV. The total number of random events was 1.57 billion, among which 1.06 billion events (67.1% of the total randoms) contained at least one photon that originated from the ¹⁷⁶Lu background radiation.

Zoom In Zoom Out Reset image size

For the anthropomorphic phantom simulation, a total of 3.17 billion coincident events were recorded (MBRD  =  20), including 35.3% trues, 15.1% scatters and 49.6% randoms. We first compared our list-mode reconstruction with different levels of data correction: (1) Trues  +  Scatters  +  Randoms without any correction (T  +  S  +  R w/o c); (2) True events only with resolution modeling (T-only w/ RM); (3) Trues  +  Scatters  +  Randoms with scatter and random correction but without resolution modeling (T  +  S  +  R w/ sc  +  rc w/o RM); and 4) Trues  +  Scatters  +  Randoms with scatter and random correction and resolution modeling (T  +  S  +  R w/ sc  +  rc  +  RM). Figure 13 and figure 14 show the reconstruction results. We can see that with scatter correction, random correction and resolution modeling (T  +  S  +  R w/sc  +  rc  +  RM), the image quality is close to the result of the trues-only reconstruction (T-only w/ RM) and more details in the lower spinal region can be seen compared with either no correction or partial correction. For quantitative comparison, figure 15 compares the CRC versus STD curves for the liver and lung lesions. The lowest curve in each group is from the reconstruction without any correction. The second lowest curve is the result with only scatter and random corrections. The two highest curves are obtained by using resolution modeling. The fact that the curve of T  +  S  +  R w/sc  +  rc  +  RM reaches almost the same CRC as that of the trues-only reconstruction indicates that the scatter and random corrections are fairly accurate.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 16 compares the reconstructed image of the EXPLORER with those from a 4-block-ring PET scanner for the imaging protocols listed in table 3. Obviously, the best image quality was achieved by the EXPLORER. The 4BR-multi produced a very noisy image. The lesion in the liver is difficult to identify. Limiting the FOV of the 4BR to one bed position reduced the noise, but the image is still noisier than that obtained by the EXPLORER.

Zoom In Zoom Out Reset image size

Figure 17 compares the CRC-STD curves for the lesion in the liver (the lung lesion is outside the 4BR single-bed FOV). The error bars were estimated from the 10 realizations of each scan protocol. It is shown that with the proper resolution modeling, all three protocols achieved a similar maximum CRC (i.e. approaching 1.0 after enough iterations), but the EXPLORER provides a 6.9-fold reduction in the standard deviation (27.9%) compared to the 4BR multi-bed whole-body imaging (192%), and 2-fold reduction (57.3%) compared to the 4BR single-bed imaging.

Zoom In Zoom Out Reset image size