- Email: [email protected]
Wobbling and LSF-based maximum likelihood expectation maximization reconstruction for wobbling PET
Radiation Physics and Chemistry 121 (2016) 1–9
Contents lists available at ScienceDirect
Radiation Physics and Chemistry journal homepage: www.elsevier.com/locate/radphyschem
Wobbling and LSF-based maximum likelihood expectation maximization reconstruction for wobbling PET Hang-Keun Kim a,b, Young-Don Son a,b, Dae-Hyuk Kwon c, Yohan Joo b, Zang-Hee Cho d,e,n a
Department of Biomedical Engineering, Gachon University, South Korea Neuroscience Research Institute, Gachon University, South Korea c Department of Medical System Engineering, Gwangju Institute of Science and Technology, South Korea d Advanced Institutes of Convergence Technology, Seoul National University, South Korea e Radiological Sciences & Psychiatry and Human Behavior, University of California, Irvine, CA 92697, USA b
H I G H L I G H T S
This paper WL-MLEM WL-MLEM WL-MLEM
proposed WL-MLEM algorithm for PET and demonstrated its performance. algorithm effectively combined wobbling and line spread function based MLEM. provided improvements in the spatial resolution and the PET image quality. can be easily extended to the other iterative reconstruction algorithms.
art ic l e i nf o
a b s t r a c t
Article history: Received 17 August 2015 Received in revised form 21 October 2015 Accepted 27 November 2015 Available online 27 November 2015
Positron emission tomography (PET) is a widely used imaging modality; however, the PET spatial resolution is not yet satisfactory for precise anatomical localization of molecular activities. Detector size is the most important factor because it determines the intrinsic resolution, which is approximately half of the detector size and determines the ultimate PET resolution. Detector size, however, cannot be made too small because both the decreased detection efficiency and the increased septal penetration effect degrade the image quality. A wobbling and line spread function (LSF)-based maximum likelihood expectation maximization (WL-MLEM) algorithm, which combined the MLEM iterative reconstruction algorithm with wobbled sampling and LSF-based deconvolution using the system matrix, was proposed for improving the spatial resolution of PET without reducing the scintillator or detector size. The new algorithm was evaluated using a simulation, and its performance was compared with that of the existing algorithms, such as conventional MLEM and LSF-based MLEM. Simulations demonstrated that the WLMLEM algorithm yielded higher spatial resolution and image quality than the existing algorithms. The WL-MLEM algorithm with wobbling PET yielded substantially improved resolution compared with conventional algorithms with stationary PET. The algorithm can be easily extended to other iterative reconstruction algorithms, such as maximum a priori (MAP) and ordered subset expectation maximization (OSEM). The WL-MLEM algorithm with wobbling PET may offer improvements in both sensitivity and resolution, the two most sought-after features in PET design. & 2015 Published by Elsevier Ltd.
Keywords: High-resolution imaging Iterative algorithms Positron emission tomography Reconstruction algorithms
1. Introduction Positron emission tomography (PET) is an imaging modality that is used widely in the clinical setting as well as in research, especially for noninvasive in vivo human studies (Cho et al., 1976; Ter-Pogossian et al., 1975). PET is the only device that can
n Corresponding author at: Advanced Institutes of Convergence Technology, Seoul National University, South Korea. E-mail addresses: [email protected], [email protected] (Z.-H. Cho).
http://dx.doi.org/10.1016/j.radphyschem.2015.11.026 0969-806X/& 2015 Published by Elsevier Ltd.
quantitatively measure molecular activity, such as the glucose metabolism, in vivo. The spatial resolution of PET, however, has been poor and not satisfactory for precise localization of much of the molecular activity that is essential for diagnosis and research, especially in brain imaging. Poor spatial resolution of PET imaging is mainly caused by a number of physical factors, such as detector size, penetration effect, positron range, and non-collinearity of annihilation photons (Derenzo et al., 1993). In PET, the effects of these physical factors are manifested as image blurring; therefore, the point spread function (PSF) of a PET system is defined by these factors. Among
2
H.-K. Kim et al. / Radiation Physics and Chemistry 121 (2016) 1–9
Fig. 1. Schematic of the three steps involved in simple two-point wobbling.
Fig. 2. Pseudo-code for the WL-MLEM algorithm showing the formula used for calculating the new pixel value x(jiter, p + 1) given the old value x(jiter, p) .
these factors, detector size is the key component that determines the intrinsic resolution, which is approximately half of the detector width. Therefore, detector size critically affects PET spatial resolution, and small-sized detectors are often preferred for highresolution systems. Reducing the detector size, however, results in decreased detection efficiency and increased septal penetration of the adjacent crystals, eventually leading to the degradation of image quality.
A number of detector developments have been proposed for high-resolution PET systems, for example, depth of interaction (DOI), which theoretically minimizes the penetration effect. However, a small detector with DOI technology is rather complex to be of practical use given the existing technology because of the increased number of channels. Under-sampling has also remained the basic problem that is responsible for the degraded spatial resolution in stationary PET (Bohm et al., 1978; Brooks et al., 1979; Herman, 1979; Huang et al., 1979; Thompson et al., 2005). Although most of the current PET systems employ interleaved sampling, which doubles the radial sampling together with half-angular sampling to meet the Nyquist sampling criterion (d/4), the situation remains far from satisfactory (Thompson et al., 2005). In the early days of PET development, various sampling techniques (Cho et al., 1983), such as halfangle rotation, wobbling, and dichotomic sampling schemes, were proposed rather than the above-described approaches to meet the requirements of the Nyquist sampling criteria. Among the various sampling schemes, wobbling sampling received great attention. However, the limited precision of the mechanical movement has been of hindrance to the practical application of the wobbling scheme. Advances in engineering and computer technology, such as precision mechanical engineering, allowed us to use wobbling in practical scanner design. In this paper, we combined the LSF deconvolution technique with wobbling by developing the wobbling and LSF-based maximum likelihood expectation maximization (WL-MLEM)
H.-K. Kim et al. / Radiation Physics and Chemistry 121 (2016) 1–9
3
Fig. 3. Translation in 2D wobbling. A wobbling PET can be easily modeled by translating the object (or the estimated image). The sinogram measured from the moved PET system (red arrows) is equivalent to the PET in which the object is moved (blue arrows). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
p (n)stationary = =
∞
∞
∞
∞
∫−∞ ∫−∞ I ( x′, y′) { δ ( x′ − nD) *LSF ( x′) } dx′dy′ ∫−∞ ∫−∞ I ( x′, y′) LSF ( x′ − nD) dx′dy′
⎡ =⎢ ⎣
⎤
∞
∫−∞ I ( x′, y′) dy′⎥⎦ *LSF ( − x′)
(2)
x ′= nD
where I (x′, y′) is the object function, x′=x cos φ +y sin φ , x′ is the radial position at a view angle ϕ . That is, this is the rotated x by the angle ϕ y′ = x sin ∅ − y cos ∅ , ⋅y′ is the same as x′ but with 90° angle and it is assumed to be the integration variable along the projection ray. That is, this is the y coordinate rotated by the angle ϕ , and * is the convolution operator. The acquired PET data are simple projection sampling data that are blurred by the LSF of the detectors. If the data are acquired in a wobbling mode at the mth position among the W wobbling points, the sampling distance can be reduced by △m, W as follows: ∞
p (n, m)wobbling =
Fig. 4. A diagram of the virtual PET that was used in this simulation study.
⎡ = ⎢ ⎣ algorithm. This combination not only satisfies the Nyquist sampling criteria but also provides an opportunity to achieve resolutions beyond the intrinsic resolution limits of PET.
∞
∫−∞ ∫−∞ I ( x′, y′) LSF ( x′ − nD − △m,W ) dx′dy′ ⎤
∞
∫−∞ I ( x′, y′) dy′⎦⎥ *LSF ( −x′)
x ′= nD +△m, W
(3)
2. Methodology and materials
Eq. (3) is the same as Eq. (2), but the additional sampling by wobbling is added; thereby, the sampling distance is reduced by △m, W , where W stands for wobbling. The blur caused by the LSF can be removed by deblurring using the estimated or measured LSF, L̃ (x):
2.1. Wobbling with LSF deconvolution
p (n, m) ⊗ L ( − x)
∞
= x = nD +Δm, W
In the general projection system, an ideal sampling function S (n) at a given view is defined by a delta function as follows:
S (n) =
∞
∞
∫−∞ ∫−∞ δ ( x′ − nD) dx′dy′
(1)
where D is the sampling interval, which is the detector width. If we assume that the 2D PET imaging system is linear shiftinvariant, the sampling function is blurred by a LSF, LSF(x), of a given detector. LSF in the PET sinogram can be regarded as the PSF in a PET image. Then, the measured sinogram, p(n), at a given view, is:
∫−∞ I (x′, y′) dy′ x′= nD +Δ
m, W
(4)
where ⊗ is the deconvolution operator. Note here that the quantity max{ △m, W } should be small enough so that the data acquired from wobbling satisfy the Nyquist sampling requirements. Fig. 1 shows the differences among p (n), p (n, m), and p (n, m) ⊗ L̃ (−x). By combining the wobbling acquisition and LSF deconvolution, it becomes possible to achieve a fine sampling of the object with sampling distance reduced by a factor of W, representing the number of wobbling stop points. This implies that we can overcome the under-sampling problem in PET systems.
4
H.-K. Kim et al. / Radiation Physics and Chemistry 121 (2016) 1–9
Fig. 5. Line spread functions (LSFs) that were used in this work. (a) The LSF of the virtual PET, which was used for the sinogram data acquisition. It had a very high sampling rate that reflected the continuity of LSFs in real PET systems. (b) Accurate (or ideal) LSF used for the ideal LSF-based image reconstruction. This was nearly identical to the LSF that was used for obtaining the measured sinogram of the virtual PET.
2.2. Wobbling and LSF-based maximum likelihood expectation maximization algorithm The WL-MLEM algorithm is, therefore, a combination of MLEM, wobbled sampling in PET, and LSF deconvolution. Before we detail the WL-MLEM, let us briefly review the MLEM equation. The update equation of common MLEM is:
xnew = x old j j
1
N
∑ aij
N ∑i = 1 aij i = 1
yi M
∑k = 1 aik x old k
(5)
where is a jth pixel value of the expected new image x. xnew j is the jth pixel value of the old expected image. xold j aij is the system matrix. yi is the ith bin value of the measured sinogram N is the number of sinogram bins M is the number of image pixels A common implementation of the MLEM algorithm uses a system matrix, which reflects the geometric factor of a PET system.
In the iterative reconstruction, LSF deconvolution is achieved by using the LSF-based fine system matrix. The finer the modeled system matrix (or the more accurate the LSF or PSF modeling), the better the spatial resolution and the signal-to-noise ratio (Jakoby et al., 2009; Olesen et al., 2009; Varrone et al., 2009). LSF-based MLEM is, therefore, an MLEM algorithm (Shepp and Vardi, 1982) with a fine system matrix that incorporates geometric factors and LSF functions. The update procedure and the pseudo-code of the WL-MLEM algorithm are given in Eq. (6) and Fig. 2, respectively. Eq. (6) and Fig. 2 show the update process and the relationship between the number of iterations and the number of wobbling points in the WL-MLEM algorithm.
⎡ 1 −1⎢ x (jiter, p + 1)=x old j Tp ⎢ N ∑ ⎣ i = 1 aij
N
∑ aij i=1
⎤ ypi ⎥ M ∑k = 1 aik Tp ⎡⎣ x (kiter, p) ⎤⎦ ⎥⎦
where iter is the current iteration number.
(6)
H.-K. Kim et al. / Radiation Physics and Chemistry 121 (2016) 1–9
5
Figs. 3 and 4). Physically, wobbled sampling corresponds to moving the center of the PET system by (x, y), and it is equivalent to translating the object by ( x, y) in the 2D case. T −p1 is, therefore, needed for updating the new estimated image. Unlike the stationary (or conventional) MLEMs, the WL-MLEM performs P (the number of wobbling points) image updates during each iteration. Therefore, the effect of a single iteration of the WLMLEM is equal to the effect of one out of p iterations of the LSFbased MLEM. Theoretically, the width of the projection ray in the LSF-based system matrix can be made as narrow as the sampling interval. Note that the sampling interval can be made even smaller than the Nyquist sampling interval by using wobbled sampling. Briefly, the LSF-based system matrix now eliminates the overlap of the line of response (LOR) functions, and the projection and back-projection using the LSFbased system matrix can reduce the width of the LOR function to the sampling interval, thus satisfying the Nyquist sampling requirement. 2.3. Materials To test the WL-MLEM algorithm, a simple 2D-based numerical simulation study was performed for evaluating the potential feasibility of wobbling in PET.
Fig. 6. Original images used in the reconstructions. The simple dot phantom with an array of holes.
x(jiter, p) is a jth pixel value of the new expected image x with iteration number iter and wobble position p. aij is the system matrix considering the LSF. ypi is the ith bin's value on the sinogram measured at the pth wobble position. N is the number of sinogram bins. M is the number of image pixels. Tp is the translational operator for pth wobble position. For wobbling PET, especially with discrete wobbling, a finite number of multiple sinogram data sets are obtained. In discrete wobbling, the number of sinograms obtained is the number of wobble positions P. In the WL-MLEM, each estimated sinogram is obtained from the translational operation and projection using the LSF-based fine system matrix. Tp is used as the translational operation for the object according to the wobble position p (see
2.3.1. Simulation parameters We designed a virtual 2D PET system to generate the sinogram data, based on the diagram in Fig. 4. The system consists of a pair of rotatable detector planes. An individual detector plane has 32 detectors with detector size of 3 mm and 100% detection efficiency, and the detectors are organized so that there are no gaps between the adjacent detectors; therefore, the detector plane is 96 mm long. Interleaving sampling was not used during the sinogram acquisition. Detector size was, therefore, equal to the radial sampling interval of the measured sinogram. The number of angular samplings (or the number of views) of the measured sinogram was 50. Angular sampling of 50 views was chosen to reduce the computation time, and this number of angular samplings was sufficient for our simulation studies. The LSFs used in this virtual PET are shown in Fig. 5. To evaluate the effect of statistical counts, Poisson noises were added to all resultant sinograms obtained from this virtual 2D PET. For simplicity, most physical processes related to detecting annihilation photons, such as attenuation, scattering, and accidental coincidence, were not considered.
Table 1 The parameters of each reconstruction. Parameters MLEM Methods Simple MLEM (MLEM)
# of Sinogramsused in recon. (or # of wobbling points) 1
System matrix
Relativetotal count ratio used for recon.
SNRin single sonogram (mean/std)
Max effectivenumber of iterations
Line-projection based (w/o LSF)
1 (low) 5 (moderate)
1
200 200
√5
10 (high)
LSF-based MLEM (LSFMLEM)
1
LSF-based
1 (low) 5 (moderate)
10
1
√5
10 (high)
Wobbling and LSF based MLEM (WL-MLEM)
5
LSF-based
1 (low) 5 (moderate)
10
1
√5
10 (high)
Wobbling and LSF based MLEM (WL-MLEM)
50
LSF-based
25 (extreme high)
10
5
200
200 200 200
40 5 40 5 40 5
20 50
6
H.-K. Kim et al. / Radiation Physics and Chemistry 121 (2016) 1–9
Fig. 7. The reconstructed images which have the least MAE. Images in the upper panel were obtained for the dot phantom at the moderate dose, and images in the lower panel were obtained for the same phantom at the low dose. The total true count used for the moderate dose was 5 times higher than that for the low dose. The results obtained by using (a) the conventional (or simple-system matrix) MLEM, (b) the LSF-MLEM, and (c) the WL-MLEM. In addition, (d) is the zoomed image (top) to the area containing holes of 1.6 mm diameter and the profile comparison (bottom) on the area of 1.6 mm diameter holes among reconstructed images. The location of cut view is displayed on the top panel of (d) (red line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
2.3.2. Digital phantoms used in the simulation For gross resolution and image reconstruction quality, a simple dot phantom of an array of small holes was used in the simulations. This phantom image consisted of 960 960 pixels (96 mm 96 mm), with image resolution of 0.1 mm/pixel. The intensity ratios for three different doses, low, moderate, and high, were applied to the phantoms (Fig. 6). The ratio was 1:5:10. In the simulation, the sinogram for the moderate-dose phantom has an SNR 5 times higher compared with that for the low-dose phantom because of Poisson noise. In addition, the extremely high dose, which had 25 times higher intensity compared with the low dose, was applied to the simple dot phantom to demonstrate the high spatial resolution capability of the WL-MLEM. Note that the intensity ratio among area A, area B, and area C was constant in all doses, and the amount of injected dose cannot change this ratio. This ratio depends on the properties, such as affinity, of the injected tracer (or radio-ligand), such as FDG, Pittsburgh compound B (PIB), and raclopride. 2.3.3. Reconstruction parameters and methods The acquired sinograms were used to reconstruct the images by using three different MLEM algorithms: the conventional MLEM, the MLEM with LSF-based system matrix (LSF-MLEM), and wobbling and LSF based MLEM. The first two MLEM methods used only a single sinogram data set that was acquired in the stationary PET mode, whereas WL-MLEM used five sinogram data sets (which were obtained at five different wobbling points). In order to set the same total counts in each algorithm, WL-MLEM used 1/5 the counts of the other algorithms' for each wobbling positions. The diameter of the wobbling trajectory was equal to the detector width, 3 mm. For every algorithm, all
voxels were set to unity in the initial iteration. The parameters that were used for the reconstruction are summarized in Table 1. It should be noted that the max effective iteration number of the algorithms was the same for all. In the WL-MLEM, multiple image updates were performed during a single iteration; therefore, the effective number of iterations should be the number of iterations multiplied by the number of wobbling points. In Table 1, therefore, the WL-MLEM with 40 iterations and five wobbling points (P¼5) corresponds to the conventional MLEM with 200 iterations. Fig. 5(b) shows the LSF kernel used in the simulation. The LSF kernel had the same shape as the ideal LSF of the virtual PET system, but with a different sampling interval. The mean absolute error (MAE) was used as a metric for the resultant reconstructed images. All reconstructed images were bilinearly interpolated before the MAE calculation in order to match the image dimension to the digital phantom (Fig. 6): M
MAE = Cdose
∑ j = 1 x recon − x original j j
M
(7)
where is a jth pixel value of the resultant reconstructed image x. x recon j is a jth pixel value of the original phantom image x. xoriginal j M is the number of image pixels. Cdose is the coefficient to normalize MAE with different doses: 1 for the low dose, 0.2 for the moderate dose, and 0.01 for the high dose. As the number of effective iteration varies in a range of 1 to 200 effective iteration, MAE were obtained for all algorithms. In order to visualize the resultant spatial resolution of each algorithm, the effective iteration number with the least MAE was chosen as the optimal images for each algorithm. The selected images were displayed and the profiles of the cross-sections were compared.
H.-K. Kim et al. / Radiation Physics and Chemistry 121 (2016) 1–9
Fig. 8. Mean absolute error (MAE) of the reconstructions. Accurate LSF was used in the reconstruction. The MAE curve in the top panel is for low dose, the middle panel for the moderate doses, and the bottom panel for the high dose. The blue solid line curve is for the conventional MLEM, the red dashed curve is for the LSFMLEM, and the red triangles are for the WL-MLEM. WL-MLEM has 5 wobble sampling positions, so that there is a single resultant image per 5 effective iterations. In the WL-MLEM (red triangle), iteration denotes the number of effective iterations (40 5), so the maximum is 200 for 5 wobbling points. The effective iteration number of the least MAE for every reconstruction: the low dose (MLEM: 14, LSF-MLEM: 13, WL-MLEM: 15), the moderate dose (MLEM:16, LSF-MLEM: 27, WL-MLEM: 30), and the high dose (MLEM:17, LSF-MLEM: 51, WL-MLEM: 35). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
3. Results PET images of an array of dot phantom were obtained by using the MLEM, LSF-MLEM, and WL-MLEM reconstruction algorithms (Fig. 7). These reconstructed images were obtained at a moderate does (upper row) and a low dose (lower row). The number of iterations was optimized in each case based on the MAE. The number of iteration which minimized MAE was selected so that the displayed resultant images were the best results for each reconstruction algorithm. Note that interleaved sampling was not used in this reconstruction. Thus, the radial sampling interval was equivalent to the detector size. As shown in Fig. 7(a), the radial sampling interval of 3 mm was too large to have good image quality in the conventional MLEM. The LSF-MLEM yielded significantly improved image quality compared with the conventional MLEM, but it could not resolve the holes that were less than 2.4 mm diameter, as seen in Fig. 7(b). Meanwhile, WLMLEM yielded the best spatial resolution, as shown in Fig. 7(c). It could partly discriminate the holes up to 1.6 mm in diameter (orange arrows). At the same dose, the WL-MLEM with five
7
wobble points significantly improved the spatial resolution and the quality images. Fig. 8 shows the MAE curves that were obtained using the MLEM, LSF-MLEM, and WL-MLEM reconstructions of the dot phantom. The accurate LSF kernel was used with LSF-MLEM and WL-MLEM. As shown in Fig. 8, the WL-MLEM showed the least MAE in all cases regardless of the dose. In particular, WL-MLEM at a moderate dose was better than LSF-MLEM at a high dose. Fig. 8 (a) shows that the accuracy of the reconstructed image improved in the WL-MLEM as the radioactive dose increased. These improvements were very limited in the MLEM and the LSF-MLEM, and it was clearly shown in the result of LSF-MLEM at high dose (Fig. 8(a)-bottom). It is because the sampling interval of LSF-MLEM was fixed as 3 mm. In contrast, image quality of the WL-MLEM improved continuously as the dose increased. Note that the conventional MLEM exhibited different MAE curves compared with other algorithms because the measured sinogram was degraded by the LSF of the PET system (Fig. 5(a)), and the bin size of the sinogram was larger than the pixel size of the image. Fig. 9 shows the spatial resolution of the WL-MLEM with higher dose activity, which simulated a PET with high sensitivity. WLMLEM with 50 sampling positions discriminated the holes of 1.2 mm diameter as shown in Fig. 9(b) and (e). Although 2.5 times increased total event count was used, the number of sampling positions increased to 50. In the view of a single sinogram, the event count decreased to 1/4 ( ¼2.5 5/50). SNR (¼ mean/std.) of a single sinogram decreased to 1/ 2 of the case of the case of high dose. As shown in Fig. 9(f), the measured single sinogram (red) had a poorer profile than did ideal sinogram (blue) because of added Poisson noise and the large detector width of 3 mm. However, WL-MLEM produced the accurate resultant image shown in Fig. 9(b) from these poor 50 sinograms. When the event count and the number of sampling positions are sufficient, the WL-MLEM can yield a high-resolution image, which has not been obtained before using any known conventional methods. A large number of (discrete) wobbling points, however, required a high enough dose (or a high system sensitivity of PET).
4. Discussion In conventional PET, the spatial resolution is determined by the intrinsic resolution and other degradation factors as discussed earlier. The PSF-MLEM algorithm(Olesen et al., 2009; Varrone et al., 2009) could remove the blurring of the PET image by deblurring using the estimated or measured PSF of the PET image. The PSF could be used as the spatial invariant or varying model. Although the PSF can be precisely measured and, thus, can eliminate all degradation factors, the spatial resolution is limited by the intrinsic resolution of the detector. The proposed WL-MLEM algorithm can substantially improve the spatial resolution of the PET image by combining the wobbling sampling strategy with the LSF deconvolution algorithm. Theoretically, we could effectively overcome the limitation of the intrinsic resolution by adding LORs and reducing sampling distance without changing the detector configuration. Required count statistics are an important factor for PET as well as the spatial resolution. A simple approach to increase the PET count statistics without changing hardware is to increase either the radiation dose or the scan time. In order to evaluate the effect of radiation dose in the WL-MLEM algorithm, the varying amounts of uptake dose were simulated using the digitized phantom. Under the same true event counts, the WL-MLEM showed improved PET spatial resolution compared with the other two reconstruction algorithms as shown in Fig. 7. In the higher dose environment as
8
H.-K. Kim et al. / Radiation Physics and Chemistry 121 (2016) 1–9
Fig. 9. Comparison of the reconstructed images of simple dot phantom at an extreme high dose. (a) The original image and (b) the result of the WL-MLEM. (c) A sinogram with ideal line projection sampling. (d) The measured sinogram from our virtual PET, with detector width of 3 mm. (e) Superimposed cut views of the two images (a, and b) and (f). Superimposed cut views of the two sinograms (c, and d), shown for comparison. The sinogram on (d) suffers from Poisson noise and the large detector width of 3 mm.
shown in Fig. 8, the SNR of the phantom was improved, and, surprisingly, the resultant PET image could resolve the holes with 1.2 mm interspace using 3 mm PET detectors. Although the spatial resolution of the conventional MLEM was the worst of all, the convergence of the algorithm was the most superior and stable in the conventional MLEM. This might have been because of the larger voxel size and thus higher SNR. Nonetheless, the minimum value in the MAE curve of the reconstructed image was the lowest in the WL-MLEM. As shown in Fig. 9, the shape of the convergence curve of WL-MLEM was nearly the same as the LSF-MLEM. Both algorithms diverged more at the lower dose as the iteration increased, which means that the characteristic of the WL-MLEM convergence might be mainly attributable to the LSF-deconvolution process and the small voxel size. In this simulation study, stopping criteria for the iteration was the least MAE to find the optimal number. MAE was calculated in all iteration and the image of the iteration which has the least MAE was selected (Figs. 7 and 8). However, it is not possible in practical use, since the true value is unknown and, thereby, MAE cannot be
calculated. MAE is an error from the true value, and it can be called the true relative error. In general, an approximate relative error from the previous iteration or a statistical hypothesis based test is used as stopping criteria (Herbert, 1990; Kontaxakis and Tzanakos, 1996; Veklerov and Llacer, 1987). These approaches can be easily implemented in the WL-MLEM, which has the similar convergence characteristic as the LSF-MLEM (Rahmim et al., 2013). Nevertheless, further study may be required for the stopping criteria of the algorithm, since the finding an optimal stopping iteration number for the iterative reconstruction is important. Due to the ill-posedness of the reconstruction problem and the inconsistency of the record from noise, there has not been yet an general solution and most users relies on the empirical approaches in clinical routines (Ben Bouallègue et al., 2013). The simplest way to improve the spatial resolution in a fixed detector system is to employ the smaller scintillators. Too small a scintillator, however, causes degradation of detection performance and increases system complexity because of the increased number of electronics channels. Thus, the width of the scintillator should
H.-K. Kim et al. / Radiation Physics and Chemistry 121 (2016) 1–9
be compromised with the system performance. However, instead of changing the detector size, we proposed a new WL-MLEM reconstruction algorithm by combining the wobbling configuration and LSF deconvolution based on the ML-EM algorithm. Our simulation confirmed that the WL-MLEM might be practically useful for improving the spatial resolution of PET without reducing the detector size and, therefore, without affecting the system sensitivity and the hardware complexity. The limitation of the proposed algorithm is that the additional hardware that drives the precise wobbling motion, is required from the conventional PET systems. Nonetheless, this newly modified algorithm suggests another recalled option to improve the system performance in addition to the recent new PET technologies, such as TOF and DOI. Another limitation of this study is that the wobbling parameters were not optimized to maximize its performance. Although Poisson noise and different doses were considered, we fixed a number of wobbling parameters, such as the number of wobbling points and wobbling diameters, in order to test the feasibility of the WL-MLEM algorithm rather than the wobbling performance. In addition, the proposed WL-MLEM can be easily extended to other reconstruction algorithms, such as maximum a posteriori (MAP) reconstruction and ordered subset EM (OSEM) algorithm.
5. Conclusion We proposed a new WL-MLEM algorithm that combined wobbling and LSF-based MLEM. This combination provided impressive improvement in the spatial resolution and the image quality with the same dose, and, thereby, it will also provide a fundamental remedy for the incurable under-sampling problems in PET. WL-MLEM reminds users of the almost forgotten technique of wobbling. In conclusion, WL-MLEM and a new PET system with precise wobbling sampling and accurate LSF modeling may be another approach for developing PET with high spatial resolution.
Acknowledgment This research was also supported by Converging Research Center Program (2013K000343) and Basic Science Research Program (NRF-2015R1C1A1A02037686) through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT & Future Planning, Korea.
References Ben Bouallègue, F., Crouzet, J.F., Mariano-Goulart, D., 2013. A heuristic statistical stopping rule for iterative reconstruction in emission tomography. Ann. Nucl. Med. 27, 84–95. http://dx.doi.org/10.1007/s12149-012-0657-5.
9
Bohm, C., Eriksson, L., Bergstrom, M., Litton, J., Sundman, R., Singh, M., 1978. A computer assisted ringdector positron camera system for reconstruction tomography of the brain. IEEE Trans. Nucl. Sci. 25, 624–637. http://dx.doi.org/ 10.1109/TNS.1978.4329384. Brooks, R.A., Sank, V.J., Talbert, A.J., Di Chiro, G., 1979. Sampling requirements and detector motion for positron emission tomography. IEEE Trans. Nucl. Sci. 26, 2760–2763. http://dx.doi.org/10.1109/TNS.1979.4330531. Cho, Z.H., Chan, J.K., Eriksson, L., 1976. Circular ring transverse axial positron camera for 3-dimensional reconstruction of radionuclides distribution. IEEE Trans. Nucl. Sci. 23, 613–622. http://dx.doi.org/10.1109/TNS.1976.4328315. Cho, Z.H., Hilal, S.K., Ra, J.B., Hong, K.S., Bigler, R.E., Yoshizumi, T., Wolf, A.P., Fowler, J.S., 1983. High-resolution circular ring positron tomograph with dichotomic sampling: dichotom-I. Phys. Med. Biol. 28, 1219. http://dx.doi.org/10.1088/ 0031-9155/28/11/003. Derenzo, S.E., Moses, W.W., Huesman, R.H., Budinger, T.F., 1993. Critical instrumentation issues for resolution smaller than 2 mm, high sensitivity brain PET. Quantif. Brain Funct. Elsevier Sci. Publ., Amst. Neth., pp. 25–37. Herbert, T.J., 1990. Statistical stopping criteria for iterative maximum likelihood reconstruction of emission images. Phys. Med. Biol. 35, 1221–1232. http://dx. doi.org/10.1088/0031-9155/35/9/003. Herman, G.T., 1979. The mathematics of wobbling a ring of positron annihilation detectors. IEEE Trans. Nucl. Sci. 26, 2756–2759. http://dx.doi.org/10.1109/ TNS.1979.4330530. Huang, S.C., Hoffman, E.J., Phelps, M.E., Kuhl, D.E., 1979. Sampling requirements of emission computed tomographic (ECT) scanners. In: 1979 The Journal of Nuclear Medicine Proceedings of the 26th Annual Meeting. Presented at the 1979 The Journal of Nuclear Medicine Proceedings of the 26th Annual Meeting, pp. 609–609. Jakoby, B.W., Bercier, Y., Watson, C.C., Bendriem, B., Townsend, D.W., 2009. Performance characteristics of a new LSO PET/CT scanner with extended axial field-of-view and PSF reconstruction. IEEE Trans. Nucl. Sci. 56, 633–639. http: //dx.doi.org/10.1109/TNS.2009.2015764. Kontaxakis, G.N., Tzanakos, G.S., 1996. Stopping criterion for the iterative EM-MLE image reconstruction for PET. Presented at the Medical Imaging 1996. Int. Soc. Opt. Photon-, 133–144. http://dx.doi.org/10.1117/12.237917. Olesen, O.V., Sibomana, M., Keller, S.H., Andersen, F., Jensen, J., Holm, S., Svarer, C., Hjgaard, L., 2009. Spatial resolution of the HRRT PET scanner using 3D OSEM PSF reconstruction. In: 2009 IEEE Nuclear Science Symposium Conference Record (NSS/MIC). Presented at the 2009 IEEE Nuclear Science Symposium Conference Record (NSS/MIC), pp. 3789–3790. doi:10.1109/NSSMIC.2009.5401892. Rahmim, A., Qi, J., Sossi, V., 2013. Resolution modeling in PET imaging: theory, practice, benefits, and pitfalls. Med. Phys. 40, 064301. http://dx.doi.org/10.1118/ 1.4800806. Shepp, L., Vardi, Y., 1982. Maximum likelihood reconstruction for emission tomography. IEEE Trans. Med. Imaging 1, 113–122. http://dx.doi.org/10.1109/ TMI.1982.4307558. Ter-Pogossian, M.M., Phelps, M.E., Hoffman, E.J., Mullani, N.A., 1975. A PositronEmission Transaxial Tomograph for Nuclear Imaging (PETT). Radiology 114, 89–98. http://dx.doi.org/10.1148/114.1.89. Thompson, C.J., James, S.S., Tomic, N., 2005. Under-sampling in PET scanners as a source of image blurring. Nucl. Instrum. Methods Phys. Res. Sect. Accel. Spectrometers Detect. Assoc. Equip. 545, 436–445. http://dx.doi.org/10.1016/j. nima.2005.01.329. Varrone, A., Sjöholm, N., Eriksson, L., Gulyás, B., Halldin, C., Farde, L., 2009. Advancement in PET quantification using 3D OP-OSEM point spread function reconstruction with the HRRT. Eur. J. Nucl. Med. Mol. Imaging 36, 1639–1650. http://dx.doi.org/10.1007/s00259-009-1156-3. Veklerov, E., Llacer, J., 1987. Stopping rule for the MLE algorithm based on statistical hypothesis testing. Med. Imaging IEEE Trans. 6, 313–319. http://dx.doi.org/ 10.1109/TMI.1987.4307849.
Contents lists available at ScienceDirect
Radiation Physics and Chemistry journal homepage: www.elsevier.com/locate/radphyschem
Wobbling and LSF-based maximum likelihood expectation maximization reconstruction for wobbling PET Hang-Keun Kim a,b, Young-Don Son a,b, Dae-Hyuk Kwon c, Yohan Joo b, Zang-Hee Cho d,e,n a
Department of Biomedical Engineering, Gachon University, South Korea Neuroscience Research Institute, Gachon University, South Korea c Department of Medical System Engineering, Gwangju Institute of Science and Technology, South Korea d Advanced Institutes of Convergence Technology, Seoul National University, South Korea e Radiological Sciences & Psychiatry and Human Behavior, University of California, Irvine, CA 92697, USA b
H I G H L I G H T S
This paper WL-MLEM WL-MLEM WL-MLEM
proposed WL-MLEM algorithm for PET and demonstrated its performance. algorithm effectively combined wobbling and line spread function based MLEM. provided improvements in the spatial resolution and the PET image quality. can be easily extended to the other iterative reconstruction algorithms.
art ic l e i nf o
a b s t r a c t
Article history: Received 17 August 2015 Received in revised form 21 October 2015 Accepted 27 November 2015 Available online 27 November 2015
Positron emission tomography (PET) is a widely used imaging modality; however, the PET spatial resolution is not yet satisfactory for precise anatomical localization of molecular activities. Detector size is the most important factor because it determines the intrinsic resolution, which is approximately half of the detector size and determines the ultimate PET resolution. Detector size, however, cannot be made too small because both the decreased detection efficiency and the increased septal penetration effect degrade the image quality. A wobbling and line spread function (LSF)-based maximum likelihood expectation maximization (WL-MLEM) algorithm, which combined the MLEM iterative reconstruction algorithm with wobbled sampling and LSF-based deconvolution using the system matrix, was proposed for improving the spatial resolution of PET without reducing the scintillator or detector size. The new algorithm was evaluated using a simulation, and its performance was compared with that of the existing algorithms, such as conventional MLEM and LSF-based MLEM. Simulations demonstrated that the WLMLEM algorithm yielded higher spatial resolution and image quality than the existing algorithms. The WL-MLEM algorithm with wobbling PET yielded substantially improved resolution compared with conventional algorithms with stationary PET. The algorithm can be easily extended to other iterative reconstruction algorithms, such as maximum a priori (MAP) and ordered subset expectation maximization (OSEM). The WL-MLEM algorithm with wobbling PET may offer improvements in both sensitivity and resolution, the two most sought-after features in PET design. & 2015 Published by Elsevier Ltd.
Keywords: High-resolution imaging Iterative algorithms Positron emission tomography Reconstruction algorithms
1. Introduction Positron emission tomography (PET) is an imaging modality that is used widely in the clinical setting as well as in research, especially for noninvasive in vivo human studies (Cho et al., 1976; Ter-Pogossian et al., 1975). PET is the only device that can
n Corresponding author at: Advanced Institutes of Convergence Technology, Seoul National University, South Korea. E-mail addresses: [email protected], [email protected] (Z.-H. Cho).
http://dx.doi.org/10.1016/j.radphyschem.2015.11.026 0969-806X/& 2015 Published by Elsevier Ltd.
quantitatively measure molecular activity, such as the glucose metabolism, in vivo. The spatial resolution of PET, however, has been poor and not satisfactory for precise localization of much of the molecular activity that is essential for diagnosis and research, especially in brain imaging. Poor spatial resolution of PET imaging is mainly caused by a number of physical factors, such as detector size, penetration effect, positron range, and non-collinearity of annihilation photons (Derenzo et al., 1993). In PET, the effects of these physical factors are manifested as image blurring; therefore, the point spread function (PSF) of a PET system is defined by these factors. Among
2
H.-K. Kim et al. / Radiation Physics and Chemistry 121 (2016) 1–9
Fig. 1. Schematic of the three steps involved in simple two-point wobbling.
Fig. 2. Pseudo-code for the WL-MLEM algorithm showing the formula used for calculating the new pixel value x(jiter, p + 1) given the old value x(jiter, p) .
these factors, detector size is the key component that determines the intrinsic resolution, which is approximately half of the detector width. Therefore, detector size critically affects PET spatial resolution, and small-sized detectors are often preferred for highresolution systems. Reducing the detector size, however, results in decreased detection efficiency and increased septal penetration of the adjacent crystals, eventually leading to the degradation of image quality.
A number of detector developments have been proposed for high-resolution PET systems, for example, depth of interaction (DOI), which theoretically minimizes the penetration effect. However, a small detector with DOI technology is rather complex to be of practical use given the existing technology because of the increased number of channels. Under-sampling has also remained the basic problem that is responsible for the degraded spatial resolution in stationary PET (Bohm et al., 1978; Brooks et al., 1979; Herman, 1979; Huang et al., 1979; Thompson et al., 2005). Although most of the current PET systems employ interleaved sampling, which doubles the radial sampling together with half-angular sampling to meet the Nyquist sampling criterion (d/4), the situation remains far from satisfactory (Thompson et al., 2005). In the early days of PET development, various sampling techniques (Cho et al., 1983), such as halfangle rotation, wobbling, and dichotomic sampling schemes, were proposed rather than the above-described approaches to meet the requirements of the Nyquist sampling criteria. Among the various sampling schemes, wobbling sampling received great attention. However, the limited precision of the mechanical movement has been of hindrance to the practical application of the wobbling scheme. Advances in engineering and computer technology, such as precision mechanical engineering, allowed us to use wobbling in practical scanner design. In this paper, we combined the LSF deconvolution technique with wobbling by developing the wobbling and LSF-based maximum likelihood expectation maximization (WL-MLEM)
H.-K. Kim et al. / Radiation Physics and Chemistry 121 (2016) 1–9
3
Fig. 3. Translation in 2D wobbling. A wobbling PET can be easily modeled by translating the object (or the estimated image). The sinogram measured from the moved PET system (red arrows) is equivalent to the PET in which the object is moved (blue arrows). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
p (n)stationary = =
∞
∞
∞
∞
∫−∞ ∫−∞ I ( x′, y′) { δ ( x′ − nD) *LSF ( x′) } dx′dy′ ∫−∞ ∫−∞ I ( x′, y′) LSF ( x′ − nD) dx′dy′
⎡ =⎢ ⎣
⎤
∞
∫−∞ I ( x′, y′) dy′⎥⎦ *LSF ( − x′)
(2)
x ′= nD
where I (x′, y′) is the object function, x′=x cos φ +y sin φ , x′ is the radial position at a view angle ϕ . That is, this is the rotated x by the angle ϕ y′ = x sin ∅ − y cos ∅ , ⋅y′ is the same as x′ but with 90° angle and it is assumed to be the integration variable along the projection ray. That is, this is the y coordinate rotated by the angle ϕ , and * is the convolution operator. The acquired PET data are simple projection sampling data that are blurred by the LSF of the detectors. If the data are acquired in a wobbling mode at the mth position among the W wobbling points, the sampling distance can be reduced by △m, W as follows: ∞
p (n, m)wobbling =
Fig. 4. A diagram of the virtual PET that was used in this simulation study.
⎡ = ⎢ ⎣ algorithm. This combination not only satisfies the Nyquist sampling criteria but also provides an opportunity to achieve resolutions beyond the intrinsic resolution limits of PET.
∞
∫−∞ ∫−∞ I ( x′, y′) LSF ( x′ − nD − △m,W ) dx′dy′ ⎤
∞
∫−∞ I ( x′, y′) dy′⎦⎥ *LSF ( −x′)
x ′= nD +△m, W
(3)
2. Methodology and materials
Eq. (3) is the same as Eq. (2), but the additional sampling by wobbling is added; thereby, the sampling distance is reduced by △m, W , where W stands for wobbling. The blur caused by the LSF can be removed by deblurring using the estimated or measured LSF, L̃ (x):
2.1. Wobbling with LSF deconvolution
p (n, m) ⊗ L ( − x)
∞
= x = nD +Δm, W
In the general projection system, an ideal sampling function S (n) at a given view is defined by a delta function as follows:
S (n) =
∞
∞
∫−∞ ∫−∞ δ ( x′ − nD) dx′dy′
(1)
where D is the sampling interval, which is the detector width. If we assume that the 2D PET imaging system is linear shiftinvariant, the sampling function is blurred by a LSF, LSF(x), of a given detector. LSF in the PET sinogram can be regarded as the PSF in a PET image. Then, the measured sinogram, p(n), at a given view, is:
∫−∞ I (x′, y′) dy′ x′= nD +Δ
m, W
(4)
where ⊗ is the deconvolution operator. Note here that the quantity max{ △m, W } should be small enough so that the data acquired from wobbling satisfy the Nyquist sampling requirements. Fig. 1 shows the differences among p (n), p (n, m), and p (n, m) ⊗ L̃ (−x). By combining the wobbling acquisition and LSF deconvolution, it becomes possible to achieve a fine sampling of the object with sampling distance reduced by a factor of W, representing the number of wobbling stop points. This implies that we can overcome the under-sampling problem in PET systems.
4
H.-K. Kim et al. / Radiation Physics and Chemistry 121 (2016) 1–9
Fig. 5. Line spread functions (LSFs) that were used in this work. (a) The LSF of the virtual PET, which was used for the sinogram data acquisition. It had a very high sampling rate that reflected the continuity of LSFs in real PET systems. (b) Accurate (or ideal) LSF used for the ideal LSF-based image reconstruction. This was nearly identical to the LSF that was used for obtaining the measured sinogram of the virtual PET.
2.2. Wobbling and LSF-based maximum likelihood expectation maximization algorithm The WL-MLEM algorithm is, therefore, a combination of MLEM, wobbled sampling in PET, and LSF deconvolution. Before we detail the WL-MLEM, let us briefly review the MLEM equation. The update equation of common MLEM is:
xnew = x old j j
1
N
∑ aij
N ∑i = 1 aij i = 1
yi M
∑k = 1 aik x old k
(5)
where is a jth pixel value of the expected new image x. xnew j is the jth pixel value of the old expected image. xold j aij is the system matrix. yi is the ith bin value of the measured sinogram N is the number of sinogram bins M is the number of image pixels A common implementation of the MLEM algorithm uses a system matrix, which reflects the geometric factor of a PET system.
In the iterative reconstruction, LSF deconvolution is achieved by using the LSF-based fine system matrix. The finer the modeled system matrix (or the more accurate the LSF or PSF modeling), the better the spatial resolution and the signal-to-noise ratio (Jakoby et al., 2009; Olesen et al., 2009; Varrone et al., 2009). LSF-based MLEM is, therefore, an MLEM algorithm (Shepp and Vardi, 1982) with a fine system matrix that incorporates geometric factors and LSF functions. The update procedure and the pseudo-code of the WL-MLEM algorithm are given in Eq. (6) and Fig. 2, respectively. Eq. (6) and Fig. 2 show the update process and the relationship between the number of iterations and the number of wobbling points in the WL-MLEM algorithm.
⎡ 1 −1⎢ x (jiter, p + 1)=x old j Tp ⎢ N ∑ ⎣ i = 1 aij
N
∑ aij i=1
⎤ ypi ⎥ M ∑k = 1 aik Tp ⎡⎣ x (kiter, p) ⎤⎦ ⎥⎦
where iter is the current iteration number.
(6)
H.-K. Kim et al. / Radiation Physics and Chemistry 121 (2016) 1–9
5
Figs. 3 and 4). Physically, wobbled sampling corresponds to moving the center of the PET system by (x, y), and it is equivalent to translating the object by ( x, y) in the 2D case. T −p1 is, therefore, needed for updating the new estimated image. Unlike the stationary (or conventional) MLEMs, the WL-MLEM performs P (the number of wobbling points) image updates during each iteration. Therefore, the effect of a single iteration of the WLMLEM is equal to the effect of one out of p iterations of the LSFbased MLEM. Theoretically, the width of the projection ray in the LSF-based system matrix can be made as narrow as the sampling interval. Note that the sampling interval can be made even smaller than the Nyquist sampling interval by using wobbled sampling. Briefly, the LSF-based system matrix now eliminates the overlap of the line of response (LOR) functions, and the projection and back-projection using the LSFbased system matrix can reduce the width of the LOR function to the sampling interval, thus satisfying the Nyquist sampling requirement. 2.3. Materials To test the WL-MLEM algorithm, a simple 2D-based numerical simulation study was performed for evaluating the potential feasibility of wobbling in PET.
Fig. 6. Original images used in the reconstructions. The simple dot phantom with an array of holes.
x(jiter, p) is a jth pixel value of the new expected image x with iteration number iter and wobble position p. aij is the system matrix considering the LSF. ypi is the ith bin's value on the sinogram measured at the pth wobble position. N is the number of sinogram bins. M is the number of image pixels. Tp is the translational operator for pth wobble position. For wobbling PET, especially with discrete wobbling, a finite number of multiple sinogram data sets are obtained. In discrete wobbling, the number of sinograms obtained is the number of wobble positions P. In the WL-MLEM, each estimated sinogram is obtained from the translational operation and projection using the LSF-based fine system matrix. Tp is used as the translational operation for the object according to the wobble position p (see
2.3.1. Simulation parameters We designed a virtual 2D PET system to generate the sinogram data, based on the diagram in Fig. 4. The system consists of a pair of rotatable detector planes. An individual detector plane has 32 detectors with detector size of 3 mm and 100% detection efficiency, and the detectors are organized so that there are no gaps between the adjacent detectors; therefore, the detector plane is 96 mm long. Interleaving sampling was not used during the sinogram acquisition. Detector size was, therefore, equal to the radial sampling interval of the measured sinogram. The number of angular samplings (or the number of views) of the measured sinogram was 50. Angular sampling of 50 views was chosen to reduce the computation time, and this number of angular samplings was sufficient for our simulation studies. The LSFs used in this virtual PET are shown in Fig. 5. To evaluate the effect of statistical counts, Poisson noises were added to all resultant sinograms obtained from this virtual 2D PET. For simplicity, most physical processes related to detecting annihilation photons, such as attenuation, scattering, and accidental coincidence, were not considered.
Table 1 The parameters of each reconstruction. Parameters MLEM Methods Simple MLEM (MLEM)
# of Sinogramsused in recon. (or # of wobbling points) 1
System matrix
Relativetotal count ratio used for recon.
SNRin single sonogram (mean/std)
Max effectivenumber of iterations
Line-projection based (w/o LSF)
1 (low) 5 (moderate)
1
200 200
√5
10 (high)
LSF-based MLEM (LSFMLEM)
1
LSF-based
1 (low) 5 (moderate)
10
1
√5
10 (high)
Wobbling and LSF based MLEM (WL-MLEM)
5
LSF-based
1 (low) 5 (moderate)
10
1
√5
10 (high)
Wobbling and LSF based MLEM (WL-MLEM)
50
LSF-based
25 (extreme high)
10
5
200
200 200 200
40 5 40 5 40 5
20 50
6
H.-K. Kim et al. / Radiation Physics and Chemistry 121 (2016) 1–9
Fig. 7. The reconstructed images which have the least MAE. Images in the upper panel were obtained for the dot phantom at the moderate dose, and images in the lower panel were obtained for the same phantom at the low dose. The total true count used for the moderate dose was 5 times higher than that for the low dose. The results obtained by using (a) the conventional (or simple-system matrix) MLEM, (b) the LSF-MLEM, and (c) the WL-MLEM. In addition, (d) is the zoomed image (top) to the area containing holes of 1.6 mm diameter and the profile comparison (bottom) on the area of 1.6 mm diameter holes among reconstructed images. The location of cut view is displayed on the top panel of (d) (red line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
2.3.2. Digital phantoms used in the simulation For gross resolution and image reconstruction quality, a simple dot phantom of an array of small holes was used in the simulations. This phantom image consisted of 960 960 pixels (96 mm 96 mm), with image resolution of 0.1 mm/pixel. The intensity ratios for three different doses, low, moderate, and high, were applied to the phantoms (Fig. 6). The ratio was 1:5:10. In the simulation, the sinogram for the moderate-dose phantom has an SNR 5 times higher compared with that for the low-dose phantom because of Poisson noise. In addition, the extremely high dose, which had 25 times higher intensity compared with the low dose, was applied to the simple dot phantom to demonstrate the high spatial resolution capability of the WL-MLEM. Note that the intensity ratio among area A, area B, and area C was constant in all doses, and the amount of injected dose cannot change this ratio. This ratio depends on the properties, such as affinity, of the injected tracer (or radio-ligand), such as FDG, Pittsburgh compound B (PIB), and raclopride. 2.3.3. Reconstruction parameters and methods The acquired sinograms were used to reconstruct the images by using three different MLEM algorithms: the conventional MLEM, the MLEM with LSF-based system matrix (LSF-MLEM), and wobbling and LSF based MLEM. The first two MLEM methods used only a single sinogram data set that was acquired in the stationary PET mode, whereas WL-MLEM used five sinogram data sets (which were obtained at five different wobbling points). In order to set the same total counts in each algorithm, WL-MLEM used 1/5 the counts of the other algorithms' for each wobbling positions. The diameter of the wobbling trajectory was equal to the detector width, 3 mm. For every algorithm, all
voxels were set to unity in the initial iteration. The parameters that were used for the reconstruction are summarized in Table 1. It should be noted that the max effective iteration number of the algorithms was the same for all. In the WL-MLEM, multiple image updates were performed during a single iteration; therefore, the effective number of iterations should be the number of iterations multiplied by the number of wobbling points. In Table 1, therefore, the WL-MLEM with 40 iterations and five wobbling points (P¼5) corresponds to the conventional MLEM with 200 iterations. Fig. 5(b) shows the LSF kernel used in the simulation. The LSF kernel had the same shape as the ideal LSF of the virtual PET system, but with a different sampling interval. The mean absolute error (MAE) was used as a metric for the resultant reconstructed images. All reconstructed images were bilinearly interpolated before the MAE calculation in order to match the image dimension to the digital phantom (Fig. 6): M
MAE = Cdose
∑ j = 1 x recon − x original j j
M
(7)
where is a jth pixel value of the resultant reconstructed image x. x recon j is a jth pixel value of the original phantom image x. xoriginal j M is the number of image pixels. Cdose is the coefficient to normalize MAE with different doses: 1 for the low dose, 0.2 for the moderate dose, and 0.01 for the high dose. As the number of effective iteration varies in a range of 1 to 200 effective iteration, MAE were obtained for all algorithms. In order to visualize the resultant spatial resolution of each algorithm, the effective iteration number with the least MAE was chosen as the optimal images for each algorithm. The selected images were displayed and the profiles of the cross-sections were compared.
H.-K. Kim et al. / Radiation Physics and Chemistry 121 (2016) 1–9
Fig. 8. Mean absolute error (MAE) of the reconstructions. Accurate LSF was used in the reconstruction. The MAE curve in the top panel is for low dose, the middle panel for the moderate doses, and the bottom panel for the high dose. The blue solid line curve is for the conventional MLEM, the red dashed curve is for the LSFMLEM, and the red triangles are for the WL-MLEM. WL-MLEM has 5 wobble sampling positions, so that there is a single resultant image per 5 effective iterations. In the WL-MLEM (red triangle), iteration denotes the number of effective iterations (40 5), so the maximum is 200 for 5 wobbling points. The effective iteration number of the least MAE for every reconstruction: the low dose (MLEM: 14, LSF-MLEM: 13, WL-MLEM: 15), the moderate dose (MLEM:16, LSF-MLEM: 27, WL-MLEM: 30), and the high dose (MLEM:17, LSF-MLEM: 51, WL-MLEM: 35). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
3. Results PET images of an array of dot phantom were obtained by using the MLEM, LSF-MLEM, and WL-MLEM reconstruction algorithms (Fig. 7). These reconstructed images were obtained at a moderate does (upper row) and a low dose (lower row). The number of iterations was optimized in each case based on the MAE. The number of iteration which minimized MAE was selected so that the displayed resultant images were the best results for each reconstruction algorithm. Note that interleaved sampling was not used in this reconstruction. Thus, the radial sampling interval was equivalent to the detector size. As shown in Fig. 7(a), the radial sampling interval of 3 mm was too large to have good image quality in the conventional MLEM. The LSF-MLEM yielded significantly improved image quality compared with the conventional MLEM, but it could not resolve the holes that were less than 2.4 mm diameter, as seen in Fig. 7(b). Meanwhile, WLMLEM yielded the best spatial resolution, as shown in Fig. 7(c). It could partly discriminate the holes up to 1.6 mm in diameter (orange arrows). At the same dose, the WL-MLEM with five
7
wobble points significantly improved the spatial resolution and the quality images. Fig. 8 shows the MAE curves that were obtained using the MLEM, LSF-MLEM, and WL-MLEM reconstructions of the dot phantom. The accurate LSF kernel was used with LSF-MLEM and WL-MLEM. As shown in Fig. 8, the WL-MLEM showed the least MAE in all cases regardless of the dose. In particular, WL-MLEM at a moderate dose was better than LSF-MLEM at a high dose. Fig. 8 (a) shows that the accuracy of the reconstructed image improved in the WL-MLEM as the radioactive dose increased. These improvements were very limited in the MLEM and the LSF-MLEM, and it was clearly shown in the result of LSF-MLEM at high dose (Fig. 8(a)-bottom). It is because the sampling interval of LSF-MLEM was fixed as 3 mm. In contrast, image quality of the WL-MLEM improved continuously as the dose increased. Note that the conventional MLEM exhibited different MAE curves compared with other algorithms because the measured sinogram was degraded by the LSF of the PET system (Fig. 5(a)), and the bin size of the sinogram was larger than the pixel size of the image. Fig. 9 shows the spatial resolution of the WL-MLEM with higher dose activity, which simulated a PET with high sensitivity. WLMLEM with 50 sampling positions discriminated the holes of 1.2 mm diameter as shown in Fig. 9(b) and (e). Although 2.5 times increased total event count was used, the number of sampling positions increased to 50. In the view of a single sinogram, the event count decreased to 1/4 ( ¼2.5 5/50). SNR (¼ mean/std.) of a single sinogram decreased to 1/ 2 of the case of the case of high dose. As shown in Fig. 9(f), the measured single sinogram (red) had a poorer profile than did ideal sinogram (blue) because of added Poisson noise and the large detector width of 3 mm. However, WL-MLEM produced the accurate resultant image shown in Fig. 9(b) from these poor 50 sinograms. When the event count and the number of sampling positions are sufficient, the WL-MLEM can yield a high-resolution image, which has not been obtained before using any known conventional methods. A large number of (discrete) wobbling points, however, required a high enough dose (or a high system sensitivity of PET).
4. Discussion In conventional PET, the spatial resolution is determined by the intrinsic resolution and other degradation factors as discussed earlier. The PSF-MLEM algorithm(Olesen et al., 2009; Varrone et al., 2009) could remove the blurring of the PET image by deblurring using the estimated or measured PSF of the PET image. The PSF could be used as the spatial invariant or varying model. Although the PSF can be precisely measured and, thus, can eliminate all degradation factors, the spatial resolution is limited by the intrinsic resolution of the detector. The proposed WL-MLEM algorithm can substantially improve the spatial resolution of the PET image by combining the wobbling sampling strategy with the LSF deconvolution algorithm. Theoretically, we could effectively overcome the limitation of the intrinsic resolution by adding LORs and reducing sampling distance without changing the detector configuration. Required count statistics are an important factor for PET as well as the spatial resolution. A simple approach to increase the PET count statistics without changing hardware is to increase either the radiation dose or the scan time. In order to evaluate the effect of radiation dose in the WL-MLEM algorithm, the varying amounts of uptake dose were simulated using the digitized phantom. Under the same true event counts, the WL-MLEM showed improved PET spatial resolution compared with the other two reconstruction algorithms as shown in Fig. 7. In the higher dose environment as
8
H.-K. Kim et al. / Radiation Physics and Chemistry 121 (2016) 1–9
Fig. 9. Comparison of the reconstructed images of simple dot phantom at an extreme high dose. (a) The original image and (b) the result of the WL-MLEM. (c) A sinogram with ideal line projection sampling. (d) The measured sinogram from our virtual PET, with detector width of 3 mm. (e) Superimposed cut views of the two images (a, and b) and (f). Superimposed cut views of the two sinograms (c, and d), shown for comparison. The sinogram on (d) suffers from Poisson noise and the large detector width of 3 mm.
shown in Fig. 8, the SNR of the phantom was improved, and, surprisingly, the resultant PET image could resolve the holes with 1.2 mm interspace using 3 mm PET detectors. Although the spatial resolution of the conventional MLEM was the worst of all, the convergence of the algorithm was the most superior and stable in the conventional MLEM. This might have been because of the larger voxel size and thus higher SNR. Nonetheless, the minimum value in the MAE curve of the reconstructed image was the lowest in the WL-MLEM. As shown in Fig. 9, the shape of the convergence curve of WL-MLEM was nearly the same as the LSF-MLEM. Both algorithms diverged more at the lower dose as the iteration increased, which means that the characteristic of the WL-MLEM convergence might be mainly attributable to the LSF-deconvolution process and the small voxel size. In this simulation study, stopping criteria for the iteration was the least MAE to find the optimal number. MAE was calculated in all iteration and the image of the iteration which has the least MAE was selected (Figs. 7 and 8). However, it is not possible in practical use, since the true value is unknown and, thereby, MAE cannot be
calculated. MAE is an error from the true value, and it can be called the true relative error. In general, an approximate relative error from the previous iteration or a statistical hypothesis based test is used as stopping criteria (Herbert, 1990; Kontaxakis and Tzanakos, 1996; Veklerov and Llacer, 1987). These approaches can be easily implemented in the WL-MLEM, which has the similar convergence characteristic as the LSF-MLEM (Rahmim et al., 2013). Nevertheless, further study may be required for the stopping criteria of the algorithm, since the finding an optimal stopping iteration number for the iterative reconstruction is important. Due to the ill-posedness of the reconstruction problem and the inconsistency of the record from noise, there has not been yet an general solution and most users relies on the empirical approaches in clinical routines (Ben Bouallègue et al., 2013). The simplest way to improve the spatial resolution in a fixed detector system is to employ the smaller scintillators. Too small a scintillator, however, causes degradation of detection performance and increases system complexity because of the increased number of electronics channels. Thus, the width of the scintillator should
H.-K. Kim et al. / Radiation Physics and Chemistry 121 (2016) 1–9
be compromised with the system performance. However, instead of changing the detector size, we proposed a new WL-MLEM reconstruction algorithm by combining the wobbling configuration and LSF deconvolution based on the ML-EM algorithm. Our simulation confirmed that the WL-MLEM might be practically useful for improving the spatial resolution of PET without reducing the detector size and, therefore, without affecting the system sensitivity and the hardware complexity. The limitation of the proposed algorithm is that the additional hardware that drives the precise wobbling motion, is required from the conventional PET systems. Nonetheless, this newly modified algorithm suggests another recalled option to improve the system performance in addition to the recent new PET technologies, such as TOF and DOI. Another limitation of this study is that the wobbling parameters were not optimized to maximize its performance. Although Poisson noise and different doses were considered, we fixed a number of wobbling parameters, such as the number of wobbling points and wobbling diameters, in order to test the feasibility of the WL-MLEM algorithm rather than the wobbling performance. In addition, the proposed WL-MLEM can be easily extended to other reconstruction algorithms, such as maximum a posteriori (MAP) reconstruction and ordered subset EM (OSEM) algorithm.
5. Conclusion We proposed a new WL-MLEM algorithm that combined wobbling and LSF-based MLEM. This combination provided impressive improvement in the spatial resolution and the image quality with the same dose, and, thereby, it will also provide a fundamental remedy for the incurable under-sampling problems in PET. WL-MLEM reminds users of the almost forgotten technique of wobbling. In conclusion, WL-MLEM and a new PET system with precise wobbling sampling and accurate LSF modeling may be another approach for developing PET with high spatial resolution.
Acknowledgment This research was also supported by Converging Research Center Program (2013K000343) and Basic Science Research Program (NRF-2015R1C1A1A02037686) through the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT & Future Planning, Korea.
References Ben Bouallègue, F., Crouzet, J.F., Mariano-Goulart, D., 2013. A heuristic statistical stopping rule for iterative reconstruction in emission tomography. Ann. Nucl. Med. 27, 84–95. http://dx.doi.org/10.1007/s12149-012-0657-5.
9
Bohm, C., Eriksson, L., Bergstrom, M., Litton, J., Sundman, R., Singh, M., 1978. A computer assisted ringdector positron camera system for reconstruction tomography of the brain. IEEE Trans. Nucl. Sci. 25, 624–637. http://dx.doi.org/ 10.1109/TNS.1978.4329384. Brooks, R.A., Sank, V.J., Talbert, A.J., Di Chiro, G., 1979. Sampling requirements and detector motion for positron emission tomography. IEEE Trans. Nucl. Sci. 26, 2760–2763. http://dx.doi.org/10.1109/TNS.1979.4330531. Cho, Z.H., Chan, J.K., Eriksson, L., 1976. Circular ring transverse axial positron camera for 3-dimensional reconstruction of radionuclides distribution. IEEE Trans. Nucl. Sci. 23, 613–622. http://dx.doi.org/10.1109/TNS.1976.4328315. Cho, Z.H., Hilal, S.K., Ra, J.B., Hong, K.S., Bigler, R.E., Yoshizumi, T., Wolf, A.P., Fowler, J.S., 1983. High-resolution circular ring positron tomograph with dichotomic sampling: dichotom-I. Phys. Med. Biol. 28, 1219. http://dx.doi.org/10.1088/ 0031-9155/28/11/003. Derenzo, S.E., Moses, W.W., Huesman, R.H., Budinger, T.F., 1993. Critical instrumentation issues for resolution smaller than 2 mm, high sensitivity brain PET. Quantif. Brain Funct. Elsevier Sci. Publ., Amst. Neth., pp. 25–37. Herbert, T.J., 1990. Statistical stopping criteria for iterative maximum likelihood reconstruction of emission images. Phys. Med. Biol. 35, 1221–1232. http://dx. doi.org/10.1088/0031-9155/35/9/003. Herman, G.T., 1979. The mathematics of wobbling a ring of positron annihilation detectors. IEEE Trans. Nucl. Sci. 26, 2756–2759. http://dx.doi.org/10.1109/ TNS.1979.4330530. Huang, S.C., Hoffman, E.J., Phelps, M.E., Kuhl, D.E., 1979. Sampling requirements of emission computed tomographic (ECT) scanners. In: 1979 The Journal of Nuclear Medicine Proceedings of the 26th Annual Meeting. Presented at the 1979 The Journal of Nuclear Medicine Proceedings of the 26th Annual Meeting, pp. 609–609. Jakoby, B.W., Bercier, Y., Watson, C.C., Bendriem, B., Townsend, D.W., 2009. Performance characteristics of a new LSO PET/CT scanner with extended axial field-of-view and PSF reconstruction. IEEE Trans. Nucl. Sci. 56, 633–639. http: //dx.doi.org/10.1109/TNS.2009.2015764. Kontaxakis, G.N., Tzanakos, G.S., 1996. Stopping criterion for the iterative EM-MLE image reconstruction for PET. Presented at the Medical Imaging 1996. Int. Soc. Opt. Photon-, 133–144. http://dx.doi.org/10.1117/12.237917. Olesen, O.V., Sibomana, M., Keller, S.H., Andersen, F., Jensen, J., Holm, S., Svarer, C., Hjgaard, L., 2009. Spatial resolution of the HRRT PET scanner using 3D OSEM PSF reconstruction. In: 2009 IEEE Nuclear Science Symposium Conference Record (NSS/MIC). Presented at the 2009 IEEE Nuclear Science Symposium Conference Record (NSS/MIC), pp. 3789–3790. doi:10.1109/NSSMIC.2009.5401892. Rahmim, A., Qi, J., Sossi, V., 2013. Resolution modeling in PET imaging: theory, practice, benefits, and pitfalls. Med. Phys. 40, 064301. http://dx.doi.org/10.1118/ 1.4800806. Shepp, L., Vardi, Y., 1982. Maximum likelihood reconstruction for emission tomography. IEEE Trans. Med. Imaging 1, 113–122. http://dx.doi.org/10.1109/ TMI.1982.4307558. Ter-Pogossian, M.M., Phelps, M.E., Hoffman, E.J., Mullani, N.A., 1975. A PositronEmission Transaxial Tomograph for Nuclear Imaging (PETT). Radiology 114, 89–98. http://dx.doi.org/10.1148/114.1.89. Thompson, C.J., James, S.S., Tomic, N., 2005. Under-sampling in PET scanners as a source of image blurring. Nucl. Instrum. Methods Phys. Res. Sect. Accel. Spectrometers Detect. Assoc. Equip. 545, 436–445. http://dx.doi.org/10.1016/j. nima.2005.01.329. Varrone, A., Sjöholm, N., Eriksson, L., Gulyás, B., Halldin, C., Farde, L., 2009. Advancement in PET quantification using 3D OP-OSEM point spread function reconstruction with the HRRT. Eur. J. Nucl. Med. Mol. Imaging 36, 1639–1650. http://dx.doi.org/10.1007/s00259-009-1156-3. Veklerov, E., Llacer, J., 1987. Stopping rule for the MLE algorithm based on statistical hypothesis testing. Med. Imaging IEEE Trans. 6, 313–319. http://dx.doi.org/ 10.1109/TMI.1987.4307849.