G Model
ARTICLE IN PRESS
CMIG-1323; No. of Pages 8
Computerized Medical Imaging and Graphics xxx (2014) xxx–xxx
Contents lists available at ScienceDirect
Computerized Medical Imaging and Graphics journal homepage: www.elsevier.com/locate/compmedimag
Acceleration of MAP-EM algorithm via over-relaxation Yu-Jung Tsai a,1 , Hsuan-Ming Huang b,1 , Yu-Hua Dean Fang c,2 , Shi-Ing Chang a,1 , Ing-Tsung Hsiao a,b,d,∗ a
Department of Medical Imaging and Radiological Sciences, College of Medicine, Chang Gung University, Taoyuan, Taiwan Medical Physics Research Center, Institute of Radiological Research, Chang Gung University and Chang Gung Memorial Hospital, Taoyuan, Taiwan Department of Electrical Engineering, Chang Gung University, Taoyuan, Taiwan d Molecular Imaging Center and Department of Nuclear Medicine, Chang Gung Memorial Hospital, Taoyuan, Taiwan b c
a r t i c l e
i n f o
Article history: Received 2 May 2013 Received in revised form 10 October 2014 Accepted 3 November 2014 Keywords: Faster tomographic reconstruction MAP-EM algorithm PET reconstruction SPECT reconstruction
a b s t r a c t To improve the convergence rate of the effective maximum a posteriori expectation-maximization (MAPEM) algorithm in tomographic reconstructions, this study proposes a modified MAP-EM which uses an over-relaxation factor to accelerate image reconstruction. The proposed method, called MAP-AEM, is evaluated and compared with the results for MAP-EM and for an ordered-subset algorithm, in terms of the convergence rate and noise properties. The results show that the proposed method converges numerically much faster than MAP-EM and with a speed that is comparable to that for an ordered-subset type method. The proposed method is effective in accelerating MAP-EM tomographic reconstruction. © 2014 Elsevier Ltd. All rights reserved.
1. Introduction In emission tomography, such as positron emission tomography (PET) and single-photon emission computed tomography (SPECT), statistical reconstructions outperform FBP (filtered backprojection) [1], because of their ability to model noise, system geometry and imaging physics, in addition to incorporating prior information that is related to the object. One of the most widely used statistical reconstruction methods is the ML-EM (maximum likelihood expectation-maximization) algorithm [2]. The ML-EM algorithm is hampered by its slow convergence rate and has difficulty with the ill-condition problem [3]. To improve its speed, the OS-EM (ordered-subset expectationmaximization) algorithm was proposed [4]. By using only one subset of the projection data per sub-iteration, the algorithm is an
∗ Corresponding author at: Department of Medical Imaging and Radiological Sciences, College of Medicine, Chang Gung University, 259 Wen-Hwa 1st Road, Kwei-Shan, Taoyuan, Taiwan. Tel.: +886 3 2118800 5399; fax: +886 3 2110052. E-mail addresses: [email protected] (Y.-J. Tsai), [email protected] (H.-M. Huang), [email protected] (Y.-H.D. Fang), [email protected] (S.-I. Chang), [email protected] (I.-T. Hsiao). 1 Address: Department of Medical Imaging and Radiological Sciences, Chang Gung University, 259 Wen-Hwa 1st Road, Kwei-Shan, Taoyuan, Taiwan. Tel.: +886 3 2118800 5389; fax: +886 3 2110052. 2 Address: Molecular Imaging Center and Nuclear Medicine, Chang Gung University and Memorial Hospital, No. 5 Fu-Shin Street, Kweishan, Taoyuan County, Taiwan. Tel.: +886 3 3281200 2625.
order of magnitude faster than the ML-EM. While OS-EM is fast, parallelizable and it preserves positivity, it does not converge to the ML solution [4]. In particular, increasing the number of subsets leads to increased noise and subset-related artifacts, which degrade lesion-detection performance [5]. Alternative methods to ensure convergence have been proposed, including RAMLA (row-action maximum likelihood algorithm) [6] and OS-separable paraboloidal surrogate [7]. RAMLA uses a manually selected relaxation parameter to preserve the property of convergence, while at the same time accelerating the iterative processing. Another way to solve the convergence problem in OSEM is the so-called complete-data ordered-subset EM (COSEM) algorithm [8], which uses all of the projection data for each pixel update, as opposed to only a subset of the projection data, which is the case for the OSEM update scheme. No extra parameter is necessary for COSEM, but the convergence rate is not as fast as that of a conventional OSEM or an optimized RAMLA. A similar reconstruction strategy, called an incremental optimization transfer algorithm, has also been applied to transmission tomography [9]. Instead of using ordered subsets, accelerated EM methods use a larger step size than the correction factor for a ML-EM [10–14]. One approach for increasing the step is to use a certain power to enlarge the correction term, while retaining the algorithm positivity [10]. By rewriting the ML-EM in an additive form, the other strategy is to multiply the correction term by a certain factor [11–14]. Recently, the COSEM reconstruction [15] used the same acceleration scheme to compensate for the slightly slow acceleration that is a result of having to process all of the data in the normalization term of the algorithm.
http://dx.doi.org/10.1016/j.compmedimag.2014.11.004 0895-6111/© 2014 Elsevier Ltd. All rights reserved.
Please cite this article in press as: Tsai Y-J, et al. Acceleration of MAP-EM algorithm via over-relaxation. Comput Med Imaging Graph (2014), http://dx.doi.org/10.1016/j.compmedimag.2014.11.004
G Model
ARTICLE IN PRESS
CMIG-1323; No. of Pages 8
Y.-J. Tsai et al. / Computerized Medical Imaging and Graphics xxx (2014) xxx–xxx
2
For the ill-condition problem, one effective solution is to add a prior or a penalty term to the ML objective function and that becomes the MAP (maximum a posteriori) reconstruction [16]. However, incorporating a prior term into a reconstruction algorithm may introduce coupling between pixels [17]. Although the OSL (one-step-late) algorithm [18] is one way to reduce the coupling in MAP-EM reconstructions, it is suboptimal and does not converge to the MAP solution [19]. As in the ML case, the convergence rate for MAP-EM is also an issue. The BSREM algorithm proposed by De Pierro and Yamagishi [19] solves the coupling problem by using a surrogate function and accelerates the reconstruction by introducing a relaxation parameter, as in the RAMLA algorithm, for convergence. Alternatively, the MAPCOSEM method, a MAP version of the COSEM algorithm [8], can be used with a similar ordered-subset updating scheme to accelerate the MAP reconstruction. In addition to the improvement in convergence rate, edge-preserving priors were developed to preserve fine detail [20–22]. However, edge enhancement that uses edgepreserving priors is limited by several factors, such as an object’s size and its intensity. To further sharpen the boundaries between tissues, while maintaining smoothness within the tissue, the incorporation of anatomical prior information in MAP reconstruction has been widely studied [23–26]. A more complete overview of MAP-EM and its related works can be found in [1]. Although the acceleration scheme [10–14] described above has been applied to conventional ML-EM, OSEM [11] and COSEM [15], it has not been applied to MAP-EM reconstruction. In this study, in order to retain good noise properties that result from the use of an entire set of projection data [11], a new algorithm that increases the convergence rate of a conventional MAP-EM is derived. The algorithm, termed MAP-AEM, incorporates the surrogate function of the quadratic prior [19] into reconstructions, without applying the OSL approach. Although these ideas were introduced at a previous conference [27], a detailed description and more simulation results, including digital and Monte Carlo simulations, are provided here. The proposed algorithms are derived in Section 2. The methods of evaluating the performance of the proposed algorithms are described in Section 3. The results, including the empirical visual performance, a comparison of the convergence rate, the reconstruction accuracy and the noise properties of the proposed methods are shown and compared with those for a conventional MAP-EM in Section 4. A comparison of the convergence rate for these methods and an OS-type reconstruction (MAP-COSEM) is also provided, to show the effectiveness of the proposed algorithms. The findings of this study are discussed in Section 5 and the conclusion follows. 2. Theory
˚P (f ) = −ˇ
the Poisson log-likelihood objective function, ˚L (g|f), and a prior ˚P (f), as: fˆ = argmax[˚L (g|f ) + ˚P (f )],
(1)
f ≥0
i
j
Hij fj +
i
gi log
j
Hij fj
(2)
(3)
where ˇ > 0 represents the smoothing parameter. Herein, N(j) is the set of neighborhood systems of pixel j and ωjj indicates the weight between pixel j and its neighboring pixel, j . Usually, √ the weight, ωjj , is 1 for horizontal and vertical neighbors and 1/ 2 for diagonal neighbors. Only 4 nearest-neighbors are used in this study. The logposterior values at the kth iteration that are used to evaluate the convergence rate are the summation of the Poisson log-likelihood objective function, ˚L (g|f), and a quadratic prior, ˚P (f), where fj is (k) replaced by the current estimate, fˆ . j
A closed-form formula is difficult to obtain because there is coupling of f in the logarithmic term of the likelihood function (2). To solve this problem, the EM algorithm was proposed [2]. It has two steps: expectation (E) and maximization (M). By substituting a surrogate function for the quadratic prior [19] and using the EM algorithm, the MAP-EM reconstruction update equations at the kth iteration can be summarized as the following two steps, where jj ≡ ωjj + ωj j and Dj ≡ i Hij [28]: (k+1)
Cij
(k) = fˆj
Hij gi
j
−b + (k+1) fˆj =
∀i, j
(4)
b2 − 4ac =d 2a
a = 4ˇ
where −
(k) Hij fˆj
j
∀j
b = −2ˇ
jj ,
(k+1) C i ij
(5)
j
(k) (k) jj (fˆj + fˆj ) + Dj ,
c=
and d = (−b + b2 − 4ac)/(2a). To derive the accelerated MAP reconstruction, a power factor, h, which is similar to that in [10–14] is introduced to these MAPEM updating equations. The accelerated MAP-EM reconstruction method is then derived by multiplying the previous update by a correction factor, which is the h power of the ratio of the new MAPEM update to the previous update, as:
(k+1) f˜j
=
(k) fˆj
h
d
(6)
(k) fˆj
By rewriting Eq. (6) in an additive form and using a Taylor series expansion, the updating equation is expressed as:
h(h − 1) 2 (k+1) (k) (k) f˜j = fˆj (1 + Xj )h = fˆj 1 + hXj + Xj + · · · 2
(7)
where (k) d − fˆj
(8)
(k) fˆj
Since the zero and first orders of a Taylor series expansion have larger values than other terms, the new updating equation is derived in terms of a first order approximation of Eq. (7) as: (k+1) (k) (k) (k) f˜j = fˆj + h(d − fˆj ) = (1 − h)fˆj + hd
(9)
In fact, the new algorithm, MAP-AEM, can be seen as an accelerated scheme that uses an additive over-relaxation factor. However, an additional count normalization term is necessary, to preserve the total counts in the data [11], as:
(k+1) (k+1) fˆj = f˜j i
˚L (g|f ) = −
2
ωjj (fj − fj ) ,
j ∈ N(j)
j
Xj =
The collected data is indicated by g and measured counts at the ith detector bin, gi (i = 1,. . ., M), have a Poisson distribution. The Poisson mean is Hf, where H is the system matrix with element Hij indicating the probability that a photon is emitted from the jth pixel and is detected by the ith bin. Also, f is the intensity object vector with pixel intensity at the jth pixel, which is expressed by fj (j = 1,. . ., N). The MAP reconstruction, fˆ , is obtained by optimizing
where
A smoothing prior, ˚P (f), of a quadratic form is used as:
g i i (k+1)
H f˜ j ij j
.
(10)
Since Eq. (10) can be combined with the projection operation of the next iteration, the additional computation time is negligible, except for the last iteration.
Please cite this article in press as: Tsai Y-J, et al. Acceleration of MAP-EM algorithm via over-relaxation. Comput Med Imaging Graph (2014), http://dx.doi.org/10.1016/j.compmedimag.2014.11.004
G Model
ARTICLE IN PRESS
CMIG-1323; No. of Pages 8
Y.-J. Tsai et al. / Computerized Medical Imaging and Graphics xxx (2014) xxx–xxx
3
Fig. 1. The digital cylindrical phantom (left), the Zubal phantom (center) and the brain phantom (right) used in this study.
3. Details of the simulations 3.1. Simulation study 1: a cylindrical phantom A 128 × 128 digital cylindrical phantom, with three regions of interest (ROIs), cold (i.e. black), background and hot (i.e. white) (with activity ratios of 1:4:8), as shown in Fig. 1(left), was used for data simulation. The system matrix, H, was generated using a ray-tracing method [29]. Poisson noise and the uniform attenuation of water at 140 keV were modeled for SPECT imaging, but not the detector response. The projection data was simulated for a total of 500 K counts, with dimensions of 128 detector bins and 128 scanning angles over 360◦ . For the noise analysis, 100 realizations of noisy data were simulated.
the striatum region to the background was 4:1. All other conditions for data simulation and acquisition were the same as those for the cylindrical phantom study. 3.3. Simulation study 3: brain phantom To investigate the performance for a more realistic condition, a brain phantom, with gray matter, white matter and cerebrospinal fluid (CSF), was simulated, as shown in Fig. 1(right). To increase complexity, three separate tumor regions were added to the phantom. The resultant intensity for tumor:gray matter:white matter:CSF is 5:3:1:0. The simulation conditions were the same as those used before. 3.4. Monte Carlo simulation
3.2. Simulation study 2: Zubal phantom To further evaluate the performance of the proposed method, a digital Zubal brain phantom [30], as shown in Fig. 1(center), was also used for simulation. To simulate the conditions for neuroimaging studies [31], the ratio of the activity concentration in
A GATE Monte Carlo simulation [32] (Version 4.0.0) was performed, to determine the feasibility of the proposed method in more realistic conditions. The platform simulates most of the physical effects during practical scanning, based on the settings for user-specified geometry, material, physics and protocol. In
Fig. 2. The reconstructed images of the cylindrical phantom, using the MAP-EM (top row) with 12, 20 and 40 iterations and using the MAP-AEM (bottom row) with h = 2 with 4, 8 and 16 iterations. The smoothing parameter is 1.0 for each reconstruction. To achieve a similar visual quality, the MAP-AEM requires fewer iterations than the MAP-EM.
Please cite this article in press as: Tsai Y-J, et al. Acceleration of MAP-EM algorithm via over-relaxation. Comput Med Imaging Graph (2014), http://dx.doi.org/10.1016/j.compmedimag.2014.11.004
G Model CMIG-1323; No. of Pages 8
ARTICLE IN PRESS Y.-J. Tsai et al. / Computerized Medical Imaging and Graphics xxx (2014) xxx–xxx
4
this study, a dual-head SPECT system, which is equipped with a typical low-energy high-resolution (LEHR) parallel-hole collimator (hole size: 1.4 mm, septal thickness: 0.16 mm, and hole length: 32.8 mm), was simulated. A NaI(Tl) crystal with an active area of 128 mm × 128 mm and thickness of 10 mm was used for photon detection. A 7.536 mL cylindrical phantom with two 0.471 mL inserts, representing cold and hot regions, was placed 47.2 mm away from the detector. The phantom contained a total of 21.5625 MBq (megabecquerels) and the relative activity concentrations for the cold, hot and background regions were the same as those seen in the simulation study 1. The cameras moved around the object over 360◦ and projection data was collected by 64 angles, taking 10 s per angle. The physical effects, including the photoelectric absorption and Compton scattering, were modeled in the simulation. An attenuation map with air (cold region) and water (hot region and background) was used. The data was then sorted into a series of sinograms, with dimension of 128 bins × 64 angles. Only the central slice, with a total of 196 K collecting counts, was used for the analysis. 3.5. Reconstructions and evaluations Following each projection data simulation, the acquired sinogram was reconstructed, using conventional MAP-EM and the proposed accelerated algorithm. For all of the reconstructions, the same system matrix, H, was the same as that used to generate the projections and an initial condition of a uniform image with a value of 1.0 was used. For the digital phantoms, a uniform attenuation map that was the same as that used in the simulation was used in the system matrix, to correct attenuation. Since the main purpose of this study is to investigate the capability of the acceleration schemes that are applied to the MAP reconstruction, the smoothing parameter, ˇ, that leads to acceptable image quality, was maintained at a fixed value (=1) for all methods, for simplicity. A larger value of ˇ leads to a smoother image and more edge information is lost, but if ˇ is too small, the reconstructed image is very noisy. Taskbased selection of the parameter is possible, but complex [33]. An empirical and easy way is to select the smoothing parameter, based on a visual investigation or the signal-to-noise ratio within some target region. However, further research is necessary to optimize the value of ˇ in an effective and efficient way. In addition, MAPAEM with various over-relaxation factors (i.e. h) was performed to investigate the increase in the speed of convergence. The performance of the proposed methods was evaluated in terms of empirical visual comparison, convergence rate, regional relative error (RE {ROI}) and noise properties, using digital phantoms. For the GATE simulation, only the visual and convergent rates are compared for the proposed and the conventional MAP-EM methods. The convergence performance is determined using the plots of log-posterior values vs. iteration numbers. Convergence is achieved when there are stable log-posterior values. All reconstructions were performed with 2, 4, 8, 16, 32 and 64 iterations. Although theoretical proof of convergence was not available, the convergence property of the proposed method was studied numerically. To investigate whether MAP-AEM converges to the same value as a convergent algorithm, such as MAP-EM, the log-posterior value of MAP-AEM at the 5000th iteration was calculated and compared with that of MAP-EM. RE {ROI} for each region was used to quantify the reconstruction accuracy of the proposed methods and the conventional MAP-EM, as: RE{ROI} ≡
1 NROI
bj
j ∈ ROI
(11)
Fig. 3. The log-posterior values for the cylindrical phantom reconstructions (MAPEM, MAP-AEM with h = 2 and for the MAP-COSEM with subset numbers of 8 and 32) under different conditions, as a function of the number of iterations, on logarithmic axes.
M 1 where bj = M |fˆ m − fj |/fj . NROI is the number of pixels within m=1 j the ROI and m is the noise trial index, for {m = 1,. . ., M}. RE {ROI} was plotted against the number of iterations and it was averaged across the 100 noise trials. A standard deviation (STD) analysis over 100 noise trials for the digital simulations was conducted to evaluate the noise properties. Similarly to the RE{ROI} evaluation, all of the results were compared with those for MAP-EM reconstruction. The STD within each ROI for each reconstruction method was computed using the following equation:
STD ≡
1 NROI
j2 ,
(12)
j ∈ ROI
M ˆ m ˆ¯ 2 M ˆ m ¯ 1 1 f . (f − f j ) and fˆ j = M where j2 = M m=1 j m=1 j The convergence rate for the proposed methods is also compared with that for the convergent OS-type MAP reconstruction and that for the MAP-COSEM with different subsets. Since using more than 32 ordered subsets did not result in significant acceleration in the preliminary experiments [15], only the ordered subsets of 8 and 32 were used. 4. Results Representative reconstructions of the cylindrical phantom for MAP-EM and MAP-AEM, with h = 2, are shown in Fig. 2, for ˇ = 1. Compared with MAP-EM, the proposed method requires a smaller number of iterations to achieve similar visual quality. Fig. 3 shows the plot of log-posterior values versus iteration numbers for MAPEM, MAP-AEM (h = 2) and MAP-COSEM (subset = 8 and 32). The MAP-AEM is faster than the MAP-EM and achieves convergence rates that are similar to those for MAP-COSEM. At the 5000th iteration, the log-posterior values for MAP-EM, MAP-AEM with h = 2, MAP-COSEM with subset = 8 and MAP-COSEM with subset = 32 are 1.454489 × 106 , 1.454487 × 106 , 1.454489 × 106 and 1.454488 × 106 , respectively. All of the reconstruction methods reach almost the same value eventually. For the 2D simulations of the cylindrical phantom, the computation time per iteration for MAP-EM, MAP-AEM with h = 2, MAP-COSEM with subset = 8 and MAP-COSEM with subset = 32, for ˇ = 1.0, are 75, 78, 125 and 138 ms, respectively. The computation
Please cite this article in press as: Tsai Y-J, et al. Acceleration of MAP-EM algorithm via over-relaxation. Comput Med Imaging Graph (2014), http://dx.doi.org/10.1016/j.compmedimag.2014.11.004
G Model CMIG-1323; No. of Pages 8
ARTICLE IN PRESS Y.-J. Tsai et al. / Computerized Medical Imaging and Graphics xxx (2014) xxx–xxx
5
Fig. 4. The average RE for the hot region (left), the cold region (middle) and the background (right) for the cylindrical phantom as a function of the number of iterations. The black (circle), blue (asterisk), green (square) and red (diamond) marks respectively represent the results for the MAP-EM, the MAP-AEM with h = 2, the MAP-COSEM with subset = 8 and the MAP-COSEM with subset = 32.
time was evaluated and reported on a personal computer with a dual Intel Pentium 4 3.00 GHz CPU. The average REs for MAP-EM, MAP-AEM (h = 2) and MAP-COSEM (subset = 8 and 32) are shown in Fig. 4 as a function of the number of iterations. Compared to MAP-EM, all of the accelerated methods achieve a lower RE after few iterations. In addition, the MAP-AEM and MAP-COSEM have similar REs, but there is some discrepancy in the early iterations. Fig. 5 shows the STD analysis for 100 noise trials with 500 K counts for the three ROIs in the cylindrical phantom. All of the accelerated methods have a higher STD than the MAP-EM in the early iterations and converge to the same point (stable STD) as the number of iterations increases. The MAP-AEM with h = 2 and the MAP-COSEM with subset = 8 show a similar acceleration, in terms of STD, but they are slower than the MAP-COSEM with subset = 32 during the early iterations. To investigate the effect of different activity distributions, the same experiments were conducted using the Zubal phantom. Fig. 6 shows anecdotal reconstructions for the MAP-EM and the MAPAEM with h = 2. Similarly to Fig. 2, all of the images are visually similar, but the MAP-AEM requires fewer iterations. Although not shown here, the acceleration, in terms of convergence rate, REs and STD studies, in reconstructing the Zubal phantom using proposed method is the same as that for the cylindrical phantom. Fig. 7 shows the reconstructed results, using both MAP-EM and MAP-AEM with h = 2 for a more realistic brain phantom. The MAP-EM requires more iterations than the MAP-AEM to achieve a similar visual quality. Fig. 8 shows the reconstructed results for a GATE simulation for the MAP-EM and the MAP-AEM with h = 2. The visual similarity of the images proves the feasibility of the proposed methods under more realistic conditions. The corresponding log-posterior value for each method, plotted against the iteration number, is
shown in Fig. 9. The result is similar to those for the digital phantom studies. 5. Discussion The capability of the acceleration scheme used in the EM methods was studied and verified for the MAP reconstruction. The digital simulations show that the proposed method outperforms the conventional MAP-EM in terms of convergence, relative error and STD stability. A comparison of the convergence speed for the proposed and the convergent OS-type MAP methods (Fig. 3) shows that the speed of convergence for the proposed method is comparable or faster if a specific over-relaxation factor, h, is used. Convergence to an optimal solution ensures reliability and stability for optimization problems. Therefore, both a rapid initial convergence rate and the ability to achieve convergence are important for clinical applications [9]. Although the most common reconstruction algorithm, OSEM, is faster than MLEM, it is non-convergent. In addition, OSEM has greater variance as the subset increases. Quantitation and lesion detection are also affected by the number of reconstruction iterations [34–36]. A fast convergent algorithm is able to provide fast, accurate and stable results. Thus, our proposed method is potentially useful to generate both fast and convergent reconstruction in clinical applications. When all of the methods use 5000 iterations to reconstruct the cylindrical phantom with a smoothing parameter of 1, they all reach similar log-posterior values. The MAP-AEM gives a slightly smaller value, possibly because there is a first-order approximation. The results for log-posterior values show that the MAP-AEM performs similarly to the MAP-EM, but it has a faster reconstruction speed, although the theoretical proof of convergence for
Fig. 5. The STD curves as a function of the number of iterations for the hot region (left), the cold region (middle) and the background (right) for the cylindrical phantom. The results for the MAP-EM (black, circle), the MAP-AEM with h = 2 (blue, asterisk), the MAP-COSEM with subset = 8 (green, square) and the MAPCOSEM with subset = 32 (red, diamond) are also shown. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
Please cite this article in press as: Tsai Y-J, et al. Acceleration of MAP-EM algorithm via over-relaxation. Comput Med Imaging Graph (2014), http://dx.doi.org/10.1016/j.compmedimag.2014.11.004
G Model CMIG-1323; No. of Pages 8 6
ARTICLE IN PRESS Y.-J. Tsai et al. / Computerized Medical Imaging and Graphics xxx (2014) xxx–xxx
Fig. 6. The anecdotal reconstructions of the Zubal phantom using the MAP-EM (top row) with 12, 20 and 40 iterations and using the MAP-AEM (bottom row) with h = 2 and for 4, 8 and 16 iterations. A smoothing parameter of 1.0 is used for all of the reconstructions. The MAP-AEM requires fewer iterations than the MAP-EM to achieve similar visual quality.
MAP-AEM is unavailable. In addition to the digital simulation, the visual comparison and the comparison of the speed of convergence for the Monte Carlo simulation show that the proposed methods accelerate MAP reconstructions under more realistic conditions. In a previous study, the value of h was approximately equal to the accelerating factor of the ML-EM algorithm [11]. The MAPAEM reconstructs faster than the MAP-EM, but this is not directly
related to the value of h. Another preliminary study (not shown here) showed that the optimal h value for the MAP-AEM is less sensitive to smoothing parameters and activity distribution. However, the optimization of h is important and a thorough investigation of the optimization of h for the MAP-AEM is worthy of future study. A complete performance comparison between the proposed methods and the OS-type MAP reconstruction is the subject of ongoing study.
Fig. 7. The reconstructed images of the digital brain phantom using the MAP-EM (top row) with 11, 20 and 54 iterations and using the MAP-AEM (bottom row) with h = 2 and with 4, 8 and 16 iterations. A smoothing parameter of 1.0 is used for all of the reconstructions. The MAP-AEM requires fewer iterations than MAP-EM to achieve similar visual quality.
Please cite this article in press as: Tsai Y-J, et al. Acceleration of MAP-EM algorithm via over-relaxation. Comput Med Imaging Graph (2014), http://dx.doi.org/10.1016/j.compmedimag.2014.11.004
G Model CMIG-1323; No. of Pages 8
ARTICLE IN PRESS Y.-J. Tsai et al. / Computerized Medical Imaging and Graphics xxx (2014) xxx–xxx
7
Fig. 8. The reconstructed images of the GATE simulation for the digital cylindrical phantom, using the MAP-EM (top row) with 12, 20 and 40 iterations and using the MAP-AEM (bottom row) with h = 2 and with 4, 8 and 16 iterations. A smoothing parameter of 1.0 is used for all of the reconstructions. The MAP-AEM requires fewer iterations than MAP-EM to achieve similar visual quality.
Although the MAP-AEM exhibits a higher convergence rate, in terms of log-posterior values, in the early iterations (Fig. 3), it is not faster in terms of RE (Fig. 4) or STD (Fig. 5). In other words, an algorithm that converges faster to reach the highest log-posterior value does not necessarily ensure optimal performance, in terms of RE or STD. This inconsistency in performance may occur because both the RE and the STD are calculated using mean regional properties and the spatial distributions of the reconstructed images in the early iterations can vary for different reconstruction methods. Nevertheless, all convergent methods eventually generate the same value. Because high-dimensional data obtained from advanced systems result in increased computational load, faster reconstruction is important for practical reasons. Since using ordered subset is one effective way to accelerate the reconstruction, future study could also include the use of an over-relaxation scheme in an OS-type MAP algorithm for further speed-up. 6. Conclusion Fig. 9. The log-posterior values of the GATE simulation using different reconstruction methods and conditions as a function of iteration numbers in a logarithmic scale by log(iteration + 1). Herein, a smoothing parameter of 1.0 was applied.
The 2D reconstructions of the cylindrical phantom show that the computation time per iteration of MAP-AEM is about 4% more than that for the MAP-EM, but the MAP-AEM is approximately twice as fast as the MAP-EM. The MAP-AEM provides similar accelerating speed performance as the MAP-COSEM does, but it requires less computation time (i.e. 62% and 57% of the computation time required for the MAP-COSEM with subset = 8 and 32, respectively). Overall, the MAP-AEM is faster than the MAP-EM and the MAP-COSEM. Note that all reconstructions were implemented in a personal computer with a dual Intel Pentium 4 3.00 GHz CPU. If implemented efficiently in a parallel system with multiprocessors, the MAP-COSEM reconstruction can be even faster.
This study proposes a method to accelerate MAP-EM reconstruction that uses an over-relaxation factor. The results show that the MAP-AEM accelerates a conventional MAP-EM reconstruction, while retaining the convergent property numerically. Using an appropriate over-relaxation factor, the proposed method even delivers a convergence rate that is comparable to that for an OStype algorithm. The proposed algorithm does not require a specific number of subsets, it eliminates the effect of a varying number of subsets and it is reliable and stable in clinical applications. Acknowledgments This work was supported by Grants MOST 103-2314-B-182010-MY3 from NSC, Taiwan, and CMRPD1C0641, CMRPD1C0642, CMRPD1C0671, CMRPD1C0672 and CMRPD1C0382 from the Research Fund of Chang Gung Memorial Hospital, Taiwan.
Please cite this article in press as: Tsai Y-J, et al. Acceleration of MAP-EM algorithm via over-relaxation. Comput Med Imaging Graph (2014), http://dx.doi.org/10.1016/j.compmedimag.2014.11.004
G Model CMIG-1323; No. of Pages 8
ARTICLE IN PRESS Y.-J. Tsai et al. / Computerized Medical Imaging and Graphics xxx (2014) xxx–xxx
8
Conflict of interest statement The authors declare that they have no conflict of interest. References [1] Qi J, Leahy RM. Iterative reconstruction techniques in emission computed tomography. Phys Med Biol 2006;51:R541–78. [2] Shepp LA, Vardi Y. Maximum likelihood reconstruction for emission tomography. IEEE Trans Med Imaging 1982;1:113–22. [3] Snyder DL, Miller MI, Thomas LJ, Politte DG. Noise and edge artifacts in maximum-likelihood reconstructions for emission tomography. IEEE Trans Med Imaging 1987;6:228–38. [4] Hudson HM, Larkin RS. Accelerated image reconstruction using ordered subsets of projection data. IEEE Trans Med Imaging 1994;13:601–9. [5] Morey AM, Kadrmas DJ. Effect of varying number of OSEM subsets on PET lesion detectability. J Nucl Med Technol 2013;41:268–73. [6] Browne J, de Pierro AB. A row-action alternative to the EM algorithm for maximizing likelihood in emission tomography. IEEE Trans Med Imaging 1996;15:687–99. [7] Ahn S, Fessler JA. Globally convergent image reconstruction for emission tomography using relaxed ordered subsets algorithms. IEEE Trans Med Imaging 2003;22:613–26. [8] Hsiao I-T, Khurd P, Rangarajan A, Gindi G. An overview of fast convergent ordered-subsets reconstruction methods for emission tomography based on the incremental EM algorithm. Nucl Instrum Methods Phys Res A 2006;569:429–33. [9] Ahn S, Fessler JA, Blatt D, Hero AO. Convergent incremental optimization transfer algorithms: application to tomography. IEEE Trans Med Imaging 2006;25:283–96. [10] Tanaka E, Nohara N, Tomitani T, Yamamoto M. Utilization of non-negativity constraints in reconstruction of emission tomograms. In: Bacharach S, editor. Inf process med imaging SE-27. Netherlands: Springer; 1986. p. 379–93. [11] Hwang D, Zeng GL. Convergence study of an accelerated ML-EM algorithm using bigger step size. Phys Med Biol 2006;51:237–52. [12] Kaufman L. Implementing and accelerating the EM algorithm for positron emission tomography. IEEE Trans Med Imaging 1987;6:37–51. [13] Lewitt RM, Muehllehner G. Accelerated iterative reconstruction for positron emission tomography based on the EM algorithm for maximum likelihood estimation. IEEE Trans Med Imaging 1986;5:16–22. [14] Tanaka E. A fast reconstruction algorithm for stationary positron emission tomography based on a modified EM algorithm. IEEE Trans Med Imaging 1987;6:98–105. [15] Hsiao I-T, Huang H-M. An accelerated ordered subsets reconstruction algorithm using an accelerating power factor for emission tomography. Phys Med Biol 2010;55:599–614. [16] Ollinger JM, Fessler JA. Positron-emission tomography. IEEE Signal Process Mag 1997;14:43–55. [17] Green PJ. Bayesian reconstructions from emission tomography data using a modified EM algorithm. IEEE Trans Med Imaging 1990;9:84–93.
[18] Green PJ. On use of the EM for penalized likelihood estimation. J R Stat Soc Ser B 1990;52:443–52. [19] De Pierro AR, Beleza Yamagishi ME. Fast EM-like methods for maximum a posteriori estimates in emission tomography. IEEE Trans Med Imaging 2001;20:280–8. [20] Alenius S, Ruotsalainen U. Bayesian image reconstruction for emission tomography based on median root prior. Eur J Nucl Med 1997;24:258–65. [21] Yu DF, Fessler JA. Edge-preserving tomographic reconstruction with nonlocal regularization. IEEE Trans Med Imaging 2002;21:159–73. [22] Hsiao I-T, Rangarajan A, Gindi G. A new convex edge-preserving median prior with applications to tomography. IEEE Trans Med Imaging 2003;22:580–5. [23] Gindi G, Lee M, Rangarajan A, Zubal IG. Bayesian reconstruction of functional images using anatomical information as priors. IEEE Trans Med Imaging 1993;12:670–80. [24] Bowsher JE, Johnson VE, Turkington TG, Jaszczak RJ, Floyd CR, Coleman RE. Bayesian reconstruction and use of anatomical a priori information for emission tomography. IEEE Trans Med Imaging 1996;15:673–86. [25] Rangarajan A, Hsiao I-T, Gindi G. A Bayesian joint mixture framework for the integration of anatomical information in functional image reconstruction. J Math Imaging Vis 2000;12:199–217. [26] Tsai Y-J, Huang H-M, Chou C-Y, Wang W, Hsiao I-T. Effective anatomical priors for emission tomographic reconstruction. J Med Biol Eng 2014, http://dx.doi.org/10.5405/jmbe.1658. Available online 13 Feb 2014. [27] Tsai Y-J, Hsiao I-T. Accelerated MAP reconstructions using an accelerated factor. In: IEEE nuclear science symposium and medical imaging conference. 2010. p. 3278–81. [28] Levitan E, Herman GT. A maximum a posteriori probability expectation maximization algorithm for image reconstruction in emission tomography. IEEE Trans Med Imaging 1987;6:185–92. [29] Siddon RL. Fast calculation of the exact radiological path for a threedimensional CT array. Med Phys 1985;12:252–5. [30] Zubal IG, Harrell CR, Smith EO, Rattner Z, Gindi G, Hoffer PB. Computerized three-dimensional segmented human anatomy. Med Phys 1994;21:299–302. [31] Lin K-J, Weng Y-H, Wey S-P, Hsiao I-T, Lu C-S, Skovronsky D, et al. Whole-body biodistribution and radiation dosimetry of 18F-FP-(+)-DTBZ (18F-AV-133): a novel vesicular monoamine transporter 2 imaging agent. J Nucl Med 2010;51:1480–5. [32] Jan S, Santin G, Strul D, Staelens S, Assié K, Autret D, et al. GATE: a simulation toolkit for PET and SPECT. Phys Med Biol 2004;49:4543–61. [33] Qi J, Huesman RH. Penalized maximum-likelihood image reconstruction for lesion detection. Phys Med Biol 2006;51:4017–29. [34] Jaskowiak CJ, Bianco JA, Perlman SB, Fine JP. Influence of reconstruction iterations on 18F-FDG PET/CT standardized uptake values. J Nucl Med 2005;46:424–8. [35] Westerterp M, Pruim J, Oyen W, Hoekstra O, Paans A, Visser E, et al. Quantification of FDG PET studies using standardised uptake values in multi-centre trials: effects of image reconstruction, resolution and ROI definition parameters. Eur J Nucl Med Mol Imaging 2007;34:392–404. [36] Inoue K, Sato T, Kitamura H, Ito M, Tsunoda Y, Hirayama A, et al. Improvement of the diagnostic accuracy of lymph node metastases of colorectal cancer in 18FFDG-PET/CT by optimizing the iteration number for the image reconstruction. Ann Nucl Med 2008;22:465–73.
Please cite this article in press as: Tsai Y-J, et al. Acceleration of MAP-EM algorithm via over-relaxation. Comput Med Imaging Graph (2014), http://dx.doi.org/10.1016/j.compmedimag.2014.11.004
ARTICLE IN PRESS
CMIG-1323; No. of Pages 8
Computerized Medical Imaging and Graphics xxx (2014) xxx–xxx
Contents lists available at ScienceDirect
Computerized Medical Imaging and Graphics journal homepage: www.elsevier.com/locate/compmedimag
Acceleration of MAP-EM algorithm via over-relaxation Yu-Jung Tsai a,1 , Hsuan-Ming Huang b,1 , Yu-Hua Dean Fang c,2 , Shi-Ing Chang a,1 , Ing-Tsung Hsiao a,b,d,∗ a
Department of Medical Imaging and Radiological Sciences, College of Medicine, Chang Gung University, Taoyuan, Taiwan Medical Physics Research Center, Institute of Radiological Research, Chang Gung University and Chang Gung Memorial Hospital, Taoyuan, Taiwan Department of Electrical Engineering, Chang Gung University, Taoyuan, Taiwan d Molecular Imaging Center and Department of Nuclear Medicine, Chang Gung Memorial Hospital, Taoyuan, Taiwan b c
a r t i c l e
i n f o
Article history: Received 2 May 2013 Received in revised form 10 October 2014 Accepted 3 November 2014 Keywords: Faster tomographic reconstruction MAP-EM algorithm PET reconstruction SPECT reconstruction
a b s t r a c t To improve the convergence rate of the effective maximum a posteriori expectation-maximization (MAPEM) algorithm in tomographic reconstructions, this study proposes a modified MAP-EM which uses an over-relaxation factor to accelerate image reconstruction. The proposed method, called MAP-AEM, is evaluated and compared with the results for MAP-EM and for an ordered-subset algorithm, in terms of the convergence rate and noise properties. The results show that the proposed method converges numerically much faster than MAP-EM and with a speed that is comparable to that for an ordered-subset type method. The proposed method is effective in accelerating MAP-EM tomographic reconstruction. © 2014 Elsevier Ltd. All rights reserved.
1. Introduction In emission tomography, such as positron emission tomography (PET) and single-photon emission computed tomography (SPECT), statistical reconstructions outperform FBP (filtered backprojection) [1], because of their ability to model noise, system geometry and imaging physics, in addition to incorporating prior information that is related to the object. One of the most widely used statistical reconstruction methods is the ML-EM (maximum likelihood expectation-maximization) algorithm [2]. The ML-EM algorithm is hampered by its slow convergence rate and has difficulty with the ill-condition problem [3]. To improve its speed, the OS-EM (ordered-subset expectationmaximization) algorithm was proposed [4]. By using only one subset of the projection data per sub-iteration, the algorithm is an
∗ Corresponding author at: Department of Medical Imaging and Radiological Sciences, College of Medicine, Chang Gung University, 259 Wen-Hwa 1st Road, Kwei-Shan, Taoyuan, Taiwan. Tel.: +886 3 2118800 5399; fax: +886 3 2110052. E-mail addresses: [email protected] (Y.-J. Tsai), [email protected] (H.-M. Huang), [email protected] (Y.-H.D. Fang), [email protected] (S.-I. Chang), [email protected] (I.-T. Hsiao). 1 Address: Department of Medical Imaging and Radiological Sciences, Chang Gung University, 259 Wen-Hwa 1st Road, Kwei-Shan, Taoyuan, Taiwan. Tel.: +886 3 2118800 5389; fax: +886 3 2110052. 2 Address: Molecular Imaging Center and Nuclear Medicine, Chang Gung University and Memorial Hospital, No. 5 Fu-Shin Street, Kweishan, Taoyuan County, Taiwan. Tel.: +886 3 3281200 2625.
order of magnitude faster than the ML-EM. While OS-EM is fast, parallelizable and it preserves positivity, it does not converge to the ML solution [4]. In particular, increasing the number of subsets leads to increased noise and subset-related artifacts, which degrade lesion-detection performance [5]. Alternative methods to ensure convergence have been proposed, including RAMLA (row-action maximum likelihood algorithm) [6] and OS-separable paraboloidal surrogate [7]. RAMLA uses a manually selected relaxation parameter to preserve the property of convergence, while at the same time accelerating the iterative processing. Another way to solve the convergence problem in OSEM is the so-called complete-data ordered-subset EM (COSEM) algorithm [8], which uses all of the projection data for each pixel update, as opposed to only a subset of the projection data, which is the case for the OSEM update scheme. No extra parameter is necessary for COSEM, but the convergence rate is not as fast as that of a conventional OSEM or an optimized RAMLA. A similar reconstruction strategy, called an incremental optimization transfer algorithm, has also been applied to transmission tomography [9]. Instead of using ordered subsets, accelerated EM methods use a larger step size than the correction factor for a ML-EM [10–14]. One approach for increasing the step is to use a certain power to enlarge the correction term, while retaining the algorithm positivity [10]. By rewriting the ML-EM in an additive form, the other strategy is to multiply the correction term by a certain factor [11–14]. Recently, the COSEM reconstruction [15] used the same acceleration scheme to compensate for the slightly slow acceleration that is a result of having to process all of the data in the normalization term of the algorithm.
http://dx.doi.org/10.1016/j.compmedimag.2014.11.004 0895-6111/© 2014 Elsevier Ltd. All rights reserved.
Please cite this article in press as: Tsai Y-J, et al. Acceleration of MAP-EM algorithm via over-relaxation. Comput Med Imaging Graph (2014), http://dx.doi.org/10.1016/j.compmedimag.2014.11.004
G Model
ARTICLE IN PRESS
CMIG-1323; No. of Pages 8
Y.-J. Tsai et al. / Computerized Medical Imaging and Graphics xxx (2014) xxx–xxx
2
For the ill-condition problem, one effective solution is to add a prior or a penalty term to the ML objective function and that becomes the MAP (maximum a posteriori) reconstruction [16]. However, incorporating a prior term into a reconstruction algorithm may introduce coupling between pixels [17]. Although the OSL (one-step-late) algorithm [18] is one way to reduce the coupling in MAP-EM reconstructions, it is suboptimal and does not converge to the MAP solution [19]. As in the ML case, the convergence rate for MAP-EM is also an issue. The BSREM algorithm proposed by De Pierro and Yamagishi [19] solves the coupling problem by using a surrogate function and accelerates the reconstruction by introducing a relaxation parameter, as in the RAMLA algorithm, for convergence. Alternatively, the MAPCOSEM method, a MAP version of the COSEM algorithm [8], can be used with a similar ordered-subset updating scheme to accelerate the MAP reconstruction. In addition to the improvement in convergence rate, edge-preserving priors were developed to preserve fine detail [20–22]. However, edge enhancement that uses edgepreserving priors is limited by several factors, such as an object’s size and its intensity. To further sharpen the boundaries between tissues, while maintaining smoothness within the tissue, the incorporation of anatomical prior information in MAP reconstruction has been widely studied [23–26]. A more complete overview of MAP-EM and its related works can be found in [1]. Although the acceleration scheme [10–14] described above has been applied to conventional ML-EM, OSEM [11] and COSEM [15], it has not been applied to MAP-EM reconstruction. In this study, in order to retain good noise properties that result from the use of an entire set of projection data [11], a new algorithm that increases the convergence rate of a conventional MAP-EM is derived. The algorithm, termed MAP-AEM, incorporates the surrogate function of the quadratic prior [19] into reconstructions, without applying the OSL approach. Although these ideas were introduced at a previous conference [27], a detailed description and more simulation results, including digital and Monte Carlo simulations, are provided here. The proposed algorithms are derived in Section 2. The methods of evaluating the performance of the proposed algorithms are described in Section 3. The results, including the empirical visual performance, a comparison of the convergence rate, the reconstruction accuracy and the noise properties of the proposed methods are shown and compared with those for a conventional MAP-EM in Section 4. A comparison of the convergence rate for these methods and an OS-type reconstruction (MAP-COSEM) is also provided, to show the effectiveness of the proposed algorithms. The findings of this study are discussed in Section 5 and the conclusion follows. 2. Theory
˚P (f ) = −ˇ
the Poisson log-likelihood objective function, ˚L (g|f), and a prior ˚P (f), as: fˆ = argmax[˚L (g|f ) + ˚P (f )],
(1)
f ≥0
i
j
Hij fj +
i
gi log
j
Hij fj
(2)
(3)
where ˇ > 0 represents the smoothing parameter. Herein, N(j) is the set of neighborhood systems of pixel j and ωjj indicates the weight between pixel j and its neighboring pixel, j . Usually, √ the weight, ωjj , is 1 for horizontal and vertical neighbors and 1/ 2 for diagonal neighbors. Only 4 nearest-neighbors are used in this study. The logposterior values at the kth iteration that are used to evaluate the convergence rate are the summation of the Poisson log-likelihood objective function, ˚L (g|f), and a quadratic prior, ˚P (f), where fj is (k) replaced by the current estimate, fˆ . j
A closed-form formula is difficult to obtain because there is coupling of f in the logarithmic term of the likelihood function (2). To solve this problem, the EM algorithm was proposed [2]. It has two steps: expectation (E) and maximization (M). By substituting a surrogate function for the quadratic prior [19] and using the EM algorithm, the MAP-EM reconstruction update equations at the kth iteration can be summarized as the following two steps, where jj ≡ ωjj + ωj j and Dj ≡ i Hij [28]: (k+1)
Cij
(k) = fˆj
Hij gi
j
−b + (k+1) fˆj =
∀i, j
(4)
b2 − 4ac =d 2a
a = 4ˇ
where −
(k) Hij fˆj
j
∀j
b = −2ˇ
jj ,
(k+1) C i ij
(5)
j
(k) (k) jj (fˆj + fˆj ) + Dj ,
c=
and d = (−b + b2 − 4ac)/(2a). To derive the accelerated MAP reconstruction, a power factor, h, which is similar to that in [10–14] is introduced to these MAPEM updating equations. The accelerated MAP-EM reconstruction method is then derived by multiplying the previous update by a correction factor, which is the h power of the ratio of the new MAPEM update to the previous update, as:
(k+1) f˜j
=
(k) fˆj
h
d
(6)
(k) fˆj
By rewriting Eq. (6) in an additive form and using a Taylor series expansion, the updating equation is expressed as:
h(h − 1) 2 (k+1) (k) (k) f˜j = fˆj (1 + Xj )h = fˆj 1 + hXj + Xj + · · · 2
(7)
where (k) d − fˆj
(8)
(k) fˆj
Since the zero and first orders of a Taylor series expansion have larger values than other terms, the new updating equation is derived in terms of a first order approximation of Eq. (7) as: (k+1) (k) (k) (k) f˜j = fˆj + h(d − fˆj ) = (1 − h)fˆj + hd
(9)
In fact, the new algorithm, MAP-AEM, can be seen as an accelerated scheme that uses an additive over-relaxation factor. However, an additional count normalization term is necessary, to preserve the total counts in the data [11], as:
(k+1) (k+1) fˆj = f˜j i
˚L (g|f ) = −
2
ωjj (fj − fj ) ,
j ∈ N(j)
j
Xj =
The collected data is indicated by g and measured counts at the ith detector bin, gi (i = 1,. . ., M), have a Poisson distribution. The Poisson mean is Hf, where H is the system matrix with element Hij indicating the probability that a photon is emitted from the jth pixel and is detected by the ith bin. Also, f is the intensity object vector with pixel intensity at the jth pixel, which is expressed by fj (j = 1,. . ., N). The MAP reconstruction, fˆ , is obtained by optimizing
where
A smoothing prior, ˚P (f), of a quadratic form is used as:
g i i (k+1)
H f˜ j ij j
.
(10)
Since Eq. (10) can be combined with the projection operation of the next iteration, the additional computation time is negligible, except for the last iteration.
Please cite this article in press as: Tsai Y-J, et al. Acceleration of MAP-EM algorithm via over-relaxation. Comput Med Imaging Graph (2014), http://dx.doi.org/10.1016/j.compmedimag.2014.11.004
G Model
ARTICLE IN PRESS
CMIG-1323; No. of Pages 8
Y.-J. Tsai et al. / Computerized Medical Imaging and Graphics xxx (2014) xxx–xxx
3
Fig. 1. The digital cylindrical phantom (left), the Zubal phantom (center) and the brain phantom (right) used in this study.
3. Details of the simulations 3.1. Simulation study 1: a cylindrical phantom A 128 × 128 digital cylindrical phantom, with three regions of interest (ROIs), cold (i.e. black), background and hot (i.e. white) (with activity ratios of 1:4:8), as shown in Fig. 1(left), was used for data simulation. The system matrix, H, was generated using a ray-tracing method [29]. Poisson noise and the uniform attenuation of water at 140 keV were modeled for SPECT imaging, but not the detector response. The projection data was simulated for a total of 500 K counts, with dimensions of 128 detector bins and 128 scanning angles over 360◦ . For the noise analysis, 100 realizations of noisy data were simulated.
the striatum region to the background was 4:1. All other conditions for data simulation and acquisition were the same as those for the cylindrical phantom study. 3.3. Simulation study 3: brain phantom To investigate the performance for a more realistic condition, a brain phantom, with gray matter, white matter and cerebrospinal fluid (CSF), was simulated, as shown in Fig. 1(right). To increase complexity, three separate tumor regions were added to the phantom. The resultant intensity for tumor:gray matter:white matter:CSF is 5:3:1:0. The simulation conditions were the same as those used before. 3.4. Monte Carlo simulation
3.2. Simulation study 2: Zubal phantom To further evaluate the performance of the proposed method, a digital Zubal brain phantom [30], as shown in Fig. 1(center), was also used for simulation. To simulate the conditions for neuroimaging studies [31], the ratio of the activity concentration in
A GATE Monte Carlo simulation [32] (Version 4.0.0) was performed, to determine the feasibility of the proposed method in more realistic conditions. The platform simulates most of the physical effects during practical scanning, based on the settings for user-specified geometry, material, physics and protocol. In
Fig. 2. The reconstructed images of the cylindrical phantom, using the MAP-EM (top row) with 12, 20 and 40 iterations and using the MAP-AEM (bottom row) with h = 2 with 4, 8 and 16 iterations. The smoothing parameter is 1.0 for each reconstruction. To achieve a similar visual quality, the MAP-AEM requires fewer iterations than the MAP-EM.
Please cite this article in press as: Tsai Y-J, et al. Acceleration of MAP-EM algorithm via over-relaxation. Comput Med Imaging Graph (2014), http://dx.doi.org/10.1016/j.compmedimag.2014.11.004
G Model CMIG-1323; No. of Pages 8
ARTICLE IN PRESS Y.-J. Tsai et al. / Computerized Medical Imaging and Graphics xxx (2014) xxx–xxx
4
this study, a dual-head SPECT system, which is equipped with a typical low-energy high-resolution (LEHR) parallel-hole collimator (hole size: 1.4 mm, septal thickness: 0.16 mm, and hole length: 32.8 mm), was simulated. A NaI(Tl) crystal with an active area of 128 mm × 128 mm and thickness of 10 mm was used for photon detection. A 7.536 mL cylindrical phantom with two 0.471 mL inserts, representing cold and hot regions, was placed 47.2 mm away from the detector. The phantom contained a total of 21.5625 MBq (megabecquerels) and the relative activity concentrations for the cold, hot and background regions were the same as those seen in the simulation study 1. The cameras moved around the object over 360◦ and projection data was collected by 64 angles, taking 10 s per angle. The physical effects, including the photoelectric absorption and Compton scattering, were modeled in the simulation. An attenuation map with air (cold region) and water (hot region and background) was used. The data was then sorted into a series of sinograms, with dimension of 128 bins × 64 angles. Only the central slice, with a total of 196 K collecting counts, was used for the analysis. 3.5. Reconstructions and evaluations Following each projection data simulation, the acquired sinogram was reconstructed, using conventional MAP-EM and the proposed accelerated algorithm. For all of the reconstructions, the same system matrix, H, was the same as that used to generate the projections and an initial condition of a uniform image with a value of 1.0 was used. For the digital phantoms, a uniform attenuation map that was the same as that used in the simulation was used in the system matrix, to correct attenuation. Since the main purpose of this study is to investigate the capability of the acceleration schemes that are applied to the MAP reconstruction, the smoothing parameter, ˇ, that leads to acceptable image quality, was maintained at a fixed value (=1) for all methods, for simplicity. A larger value of ˇ leads to a smoother image and more edge information is lost, but if ˇ is too small, the reconstructed image is very noisy. Taskbased selection of the parameter is possible, but complex [33]. An empirical and easy way is to select the smoothing parameter, based on a visual investigation or the signal-to-noise ratio within some target region. However, further research is necessary to optimize the value of ˇ in an effective and efficient way. In addition, MAPAEM with various over-relaxation factors (i.e. h) was performed to investigate the increase in the speed of convergence. The performance of the proposed methods was evaluated in terms of empirical visual comparison, convergence rate, regional relative error (RE {ROI}) and noise properties, using digital phantoms. For the GATE simulation, only the visual and convergent rates are compared for the proposed and the conventional MAP-EM methods. The convergence performance is determined using the plots of log-posterior values vs. iteration numbers. Convergence is achieved when there are stable log-posterior values. All reconstructions were performed with 2, 4, 8, 16, 32 and 64 iterations. Although theoretical proof of convergence was not available, the convergence property of the proposed method was studied numerically. To investigate whether MAP-AEM converges to the same value as a convergent algorithm, such as MAP-EM, the log-posterior value of MAP-AEM at the 5000th iteration was calculated and compared with that of MAP-EM. RE {ROI} for each region was used to quantify the reconstruction accuracy of the proposed methods and the conventional MAP-EM, as: RE{ROI} ≡
1 NROI
bj
j ∈ ROI
(11)
Fig. 3. The log-posterior values for the cylindrical phantom reconstructions (MAPEM, MAP-AEM with h = 2 and for the MAP-COSEM with subset numbers of 8 and 32) under different conditions, as a function of the number of iterations, on logarithmic axes.
M 1 where bj = M |fˆ m − fj |/fj . NROI is the number of pixels within m=1 j the ROI and m is the noise trial index, for {m = 1,. . ., M}. RE {ROI} was plotted against the number of iterations and it was averaged across the 100 noise trials. A standard deviation (STD) analysis over 100 noise trials for the digital simulations was conducted to evaluate the noise properties. Similarly to the RE{ROI} evaluation, all of the results were compared with those for MAP-EM reconstruction. The STD within each ROI for each reconstruction method was computed using the following equation:
STD ≡
1 NROI
j2 ,
(12)
j ∈ ROI
M ˆ m ˆ¯ 2 M ˆ m ¯ 1 1 f . (f − f j ) and fˆ j = M where j2 = M m=1 j m=1 j The convergence rate for the proposed methods is also compared with that for the convergent OS-type MAP reconstruction and that for the MAP-COSEM with different subsets. Since using more than 32 ordered subsets did not result in significant acceleration in the preliminary experiments [15], only the ordered subsets of 8 and 32 were used. 4. Results Representative reconstructions of the cylindrical phantom for MAP-EM and MAP-AEM, with h = 2, are shown in Fig. 2, for ˇ = 1. Compared with MAP-EM, the proposed method requires a smaller number of iterations to achieve similar visual quality. Fig. 3 shows the plot of log-posterior values versus iteration numbers for MAPEM, MAP-AEM (h = 2) and MAP-COSEM (subset = 8 and 32). The MAP-AEM is faster than the MAP-EM and achieves convergence rates that are similar to those for MAP-COSEM. At the 5000th iteration, the log-posterior values for MAP-EM, MAP-AEM with h = 2, MAP-COSEM with subset = 8 and MAP-COSEM with subset = 32 are 1.454489 × 106 , 1.454487 × 106 , 1.454489 × 106 and 1.454488 × 106 , respectively. All of the reconstruction methods reach almost the same value eventually. For the 2D simulations of the cylindrical phantom, the computation time per iteration for MAP-EM, MAP-AEM with h = 2, MAP-COSEM with subset = 8 and MAP-COSEM with subset = 32, for ˇ = 1.0, are 75, 78, 125 and 138 ms, respectively. The computation
Please cite this article in press as: Tsai Y-J, et al. Acceleration of MAP-EM algorithm via over-relaxation. Comput Med Imaging Graph (2014), http://dx.doi.org/10.1016/j.compmedimag.2014.11.004
G Model CMIG-1323; No. of Pages 8
ARTICLE IN PRESS Y.-J. Tsai et al. / Computerized Medical Imaging and Graphics xxx (2014) xxx–xxx
5
Fig. 4. The average RE for the hot region (left), the cold region (middle) and the background (right) for the cylindrical phantom as a function of the number of iterations. The black (circle), blue (asterisk), green (square) and red (diamond) marks respectively represent the results for the MAP-EM, the MAP-AEM with h = 2, the MAP-COSEM with subset = 8 and the MAP-COSEM with subset = 32.
time was evaluated and reported on a personal computer with a dual Intel Pentium 4 3.00 GHz CPU. The average REs for MAP-EM, MAP-AEM (h = 2) and MAP-COSEM (subset = 8 and 32) are shown in Fig. 4 as a function of the number of iterations. Compared to MAP-EM, all of the accelerated methods achieve a lower RE after few iterations. In addition, the MAP-AEM and MAP-COSEM have similar REs, but there is some discrepancy in the early iterations. Fig. 5 shows the STD analysis for 100 noise trials with 500 K counts for the three ROIs in the cylindrical phantom. All of the accelerated methods have a higher STD than the MAP-EM in the early iterations and converge to the same point (stable STD) as the number of iterations increases. The MAP-AEM with h = 2 and the MAP-COSEM with subset = 8 show a similar acceleration, in terms of STD, but they are slower than the MAP-COSEM with subset = 32 during the early iterations. To investigate the effect of different activity distributions, the same experiments were conducted using the Zubal phantom. Fig. 6 shows anecdotal reconstructions for the MAP-EM and the MAPAEM with h = 2. Similarly to Fig. 2, all of the images are visually similar, but the MAP-AEM requires fewer iterations. Although not shown here, the acceleration, in terms of convergence rate, REs and STD studies, in reconstructing the Zubal phantom using proposed method is the same as that for the cylindrical phantom. Fig. 7 shows the reconstructed results, using both MAP-EM and MAP-AEM with h = 2 for a more realistic brain phantom. The MAP-EM requires more iterations than the MAP-AEM to achieve a similar visual quality. Fig. 8 shows the reconstructed results for a GATE simulation for the MAP-EM and the MAP-AEM with h = 2. The visual similarity of the images proves the feasibility of the proposed methods under more realistic conditions. The corresponding log-posterior value for each method, plotted against the iteration number, is
shown in Fig. 9. The result is similar to those for the digital phantom studies. 5. Discussion The capability of the acceleration scheme used in the EM methods was studied and verified for the MAP reconstruction. The digital simulations show that the proposed method outperforms the conventional MAP-EM in terms of convergence, relative error and STD stability. A comparison of the convergence speed for the proposed and the convergent OS-type MAP methods (Fig. 3) shows that the speed of convergence for the proposed method is comparable or faster if a specific over-relaxation factor, h, is used. Convergence to an optimal solution ensures reliability and stability for optimization problems. Therefore, both a rapid initial convergence rate and the ability to achieve convergence are important for clinical applications [9]. Although the most common reconstruction algorithm, OSEM, is faster than MLEM, it is non-convergent. In addition, OSEM has greater variance as the subset increases. Quantitation and lesion detection are also affected by the number of reconstruction iterations [34–36]. A fast convergent algorithm is able to provide fast, accurate and stable results. Thus, our proposed method is potentially useful to generate both fast and convergent reconstruction in clinical applications. When all of the methods use 5000 iterations to reconstruct the cylindrical phantom with a smoothing parameter of 1, they all reach similar log-posterior values. The MAP-AEM gives a slightly smaller value, possibly because there is a first-order approximation. The results for log-posterior values show that the MAP-AEM performs similarly to the MAP-EM, but it has a faster reconstruction speed, although the theoretical proof of convergence for
Fig. 5. The STD curves as a function of the number of iterations for the hot region (left), the cold region (middle) and the background (right) for the cylindrical phantom. The results for the MAP-EM (black, circle), the MAP-AEM with h = 2 (blue, asterisk), the MAP-COSEM with subset = 8 (green, square) and the MAPCOSEM with subset = 32 (red, diamond) are also shown. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.)
Please cite this article in press as: Tsai Y-J, et al. Acceleration of MAP-EM algorithm via over-relaxation. Comput Med Imaging Graph (2014), http://dx.doi.org/10.1016/j.compmedimag.2014.11.004
G Model CMIG-1323; No. of Pages 8 6
ARTICLE IN PRESS Y.-J. Tsai et al. / Computerized Medical Imaging and Graphics xxx (2014) xxx–xxx
Fig. 6. The anecdotal reconstructions of the Zubal phantom using the MAP-EM (top row) with 12, 20 and 40 iterations and using the MAP-AEM (bottom row) with h = 2 and for 4, 8 and 16 iterations. A smoothing parameter of 1.0 is used for all of the reconstructions. The MAP-AEM requires fewer iterations than the MAP-EM to achieve similar visual quality.
MAP-AEM is unavailable. In addition to the digital simulation, the visual comparison and the comparison of the speed of convergence for the Monte Carlo simulation show that the proposed methods accelerate MAP reconstructions under more realistic conditions. In a previous study, the value of h was approximately equal to the accelerating factor of the ML-EM algorithm [11]. The MAPAEM reconstructs faster than the MAP-EM, but this is not directly
related to the value of h. Another preliminary study (not shown here) showed that the optimal h value for the MAP-AEM is less sensitive to smoothing parameters and activity distribution. However, the optimization of h is important and a thorough investigation of the optimization of h for the MAP-AEM is worthy of future study. A complete performance comparison between the proposed methods and the OS-type MAP reconstruction is the subject of ongoing study.
Fig. 7. The reconstructed images of the digital brain phantom using the MAP-EM (top row) with 11, 20 and 54 iterations and using the MAP-AEM (bottom row) with h = 2 and with 4, 8 and 16 iterations. A smoothing parameter of 1.0 is used for all of the reconstructions. The MAP-AEM requires fewer iterations than MAP-EM to achieve similar visual quality.
Please cite this article in press as: Tsai Y-J, et al. Acceleration of MAP-EM algorithm via over-relaxation. Comput Med Imaging Graph (2014), http://dx.doi.org/10.1016/j.compmedimag.2014.11.004
G Model CMIG-1323; No. of Pages 8
ARTICLE IN PRESS Y.-J. Tsai et al. / Computerized Medical Imaging and Graphics xxx (2014) xxx–xxx
7
Fig. 8. The reconstructed images of the GATE simulation for the digital cylindrical phantom, using the MAP-EM (top row) with 12, 20 and 40 iterations and using the MAP-AEM (bottom row) with h = 2 and with 4, 8 and 16 iterations. A smoothing parameter of 1.0 is used for all of the reconstructions. The MAP-AEM requires fewer iterations than MAP-EM to achieve similar visual quality.
Although the MAP-AEM exhibits a higher convergence rate, in terms of log-posterior values, in the early iterations (Fig. 3), it is not faster in terms of RE (Fig. 4) or STD (Fig. 5). In other words, an algorithm that converges faster to reach the highest log-posterior value does not necessarily ensure optimal performance, in terms of RE or STD. This inconsistency in performance may occur because both the RE and the STD are calculated using mean regional properties and the spatial distributions of the reconstructed images in the early iterations can vary for different reconstruction methods. Nevertheless, all convergent methods eventually generate the same value. Because high-dimensional data obtained from advanced systems result in increased computational load, faster reconstruction is important for practical reasons. Since using ordered subset is one effective way to accelerate the reconstruction, future study could also include the use of an over-relaxation scheme in an OS-type MAP algorithm for further speed-up. 6. Conclusion Fig. 9. The log-posterior values of the GATE simulation using different reconstruction methods and conditions as a function of iteration numbers in a logarithmic scale by log(iteration + 1). Herein, a smoothing parameter of 1.0 was applied.
The 2D reconstructions of the cylindrical phantom show that the computation time per iteration of MAP-AEM is about 4% more than that for the MAP-EM, but the MAP-AEM is approximately twice as fast as the MAP-EM. The MAP-AEM provides similar accelerating speed performance as the MAP-COSEM does, but it requires less computation time (i.e. 62% and 57% of the computation time required for the MAP-COSEM with subset = 8 and 32, respectively). Overall, the MAP-AEM is faster than the MAP-EM and the MAP-COSEM. Note that all reconstructions were implemented in a personal computer with a dual Intel Pentium 4 3.00 GHz CPU. If implemented efficiently in a parallel system with multiprocessors, the MAP-COSEM reconstruction can be even faster.
This study proposes a method to accelerate MAP-EM reconstruction that uses an over-relaxation factor. The results show that the MAP-AEM accelerates a conventional MAP-EM reconstruction, while retaining the convergent property numerically. Using an appropriate over-relaxation factor, the proposed method even delivers a convergence rate that is comparable to that for an OStype algorithm. The proposed algorithm does not require a specific number of subsets, it eliminates the effect of a varying number of subsets and it is reliable and stable in clinical applications. Acknowledgments This work was supported by Grants MOST 103-2314-B-182010-MY3 from NSC, Taiwan, and CMRPD1C0641, CMRPD1C0642, CMRPD1C0671, CMRPD1C0672 and CMRPD1C0382 from the Research Fund of Chang Gung Memorial Hospital, Taiwan.
Please cite this article in press as: Tsai Y-J, et al. Acceleration of MAP-EM algorithm via over-relaxation. Comput Med Imaging Graph (2014), http://dx.doi.org/10.1016/j.compmedimag.2014.11.004
G Model CMIG-1323; No. of Pages 8
ARTICLE IN PRESS Y.-J. Tsai et al. / Computerized Medical Imaging and Graphics xxx (2014) xxx–xxx
8
Conflict of interest statement The authors declare that they have no conflict of interest. References [1] Qi J, Leahy RM. Iterative reconstruction techniques in emission computed tomography. Phys Med Biol 2006;51:R541–78. [2] Shepp LA, Vardi Y. Maximum likelihood reconstruction for emission tomography. IEEE Trans Med Imaging 1982;1:113–22. [3] Snyder DL, Miller MI, Thomas LJ, Politte DG. Noise and edge artifacts in maximum-likelihood reconstructions for emission tomography. IEEE Trans Med Imaging 1987;6:228–38. [4] Hudson HM, Larkin RS. Accelerated image reconstruction using ordered subsets of projection data. IEEE Trans Med Imaging 1994;13:601–9. [5] Morey AM, Kadrmas DJ. Effect of varying number of OSEM subsets on PET lesion detectability. J Nucl Med Technol 2013;41:268–73. [6] Browne J, de Pierro AB. A row-action alternative to the EM algorithm for maximizing likelihood in emission tomography. IEEE Trans Med Imaging 1996;15:687–99. [7] Ahn S, Fessler JA. Globally convergent image reconstruction for emission tomography using relaxed ordered subsets algorithms. IEEE Trans Med Imaging 2003;22:613–26. [8] Hsiao I-T, Khurd P, Rangarajan A, Gindi G. An overview of fast convergent ordered-subsets reconstruction methods for emission tomography based on the incremental EM algorithm. Nucl Instrum Methods Phys Res A 2006;569:429–33. [9] Ahn S, Fessler JA, Blatt D, Hero AO. Convergent incremental optimization transfer algorithms: application to tomography. IEEE Trans Med Imaging 2006;25:283–96. [10] Tanaka E, Nohara N, Tomitani T, Yamamoto M. Utilization of non-negativity constraints in reconstruction of emission tomograms. In: Bacharach S, editor. Inf process med imaging SE-27. Netherlands: Springer; 1986. p. 379–93. [11] Hwang D, Zeng GL. Convergence study of an accelerated ML-EM algorithm using bigger step size. Phys Med Biol 2006;51:237–52. [12] Kaufman L. Implementing and accelerating the EM algorithm for positron emission tomography. IEEE Trans Med Imaging 1987;6:37–51. [13] Lewitt RM, Muehllehner G. Accelerated iterative reconstruction for positron emission tomography based on the EM algorithm for maximum likelihood estimation. IEEE Trans Med Imaging 1986;5:16–22. [14] Tanaka E. A fast reconstruction algorithm for stationary positron emission tomography based on a modified EM algorithm. IEEE Trans Med Imaging 1987;6:98–105. [15] Hsiao I-T, Huang H-M. An accelerated ordered subsets reconstruction algorithm using an accelerating power factor for emission tomography. Phys Med Biol 2010;55:599–614. [16] Ollinger JM, Fessler JA. Positron-emission tomography. IEEE Signal Process Mag 1997;14:43–55. [17] Green PJ. Bayesian reconstructions from emission tomography data using a modified EM algorithm. IEEE Trans Med Imaging 1990;9:84–93.
[18] Green PJ. On use of the EM for penalized likelihood estimation. J R Stat Soc Ser B 1990;52:443–52. [19] De Pierro AR, Beleza Yamagishi ME. Fast EM-like methods for maximum a posteriori estimates in emission tomography. IEEE Trans Med Imaging 2001;20:280–8. [20] Alenius S, Ruotsalainen U. Bayesian image reconstruction for emission tomography based on median root prior. Eur J Nucl Med 1997;24:258–65. [21] Yu DF, Fessler JA. Edge-preserving tomographic reconstruction with nonlocal regularization. IEEE Trans Med Imaging 2002;21:159–73. [22] Hsiao I-T, Rangarajan A, Gindi G. A new convex edge-preserving median prior with applications to tomography. IEEE Trans Med Imaging 2003;22:580–5. [23] Gindi G, Lee M, Rangarajan A, Zubal IG. Bayesian reconstruction of functional images using anatomical information as priors. IEEE Trans Med Imaging 1993;12:670–80. [24] Bowsher JE, Johnson VE, Turkington TG, Jaszczak RJ, Floyd CR, Coleman RE. Bayesian reconstruction and use of anatomical a priori information for emission tomography. IEEE Trans Med Imaging 1996;15:673–86. [25] Rangarajan A, Hsiao I-T, Gindi G. A Bayesian joint mixture framework for the integration of anatomical information in functional image reconstruction. J Math Imaging Vis 2000;12:199–217. [26] Tsai Y-J, Huang H-M, Chou C-Y, Wang W, Hsiao I-T. Effective anatomical priors for emission tomographic reconstruction. J Med Biol Eng 2014, http://dx.doi.org/10.5405/jmbe.1658. Available online 13 Feb 2014. [27] Tsai Y-J, Hsiao I-T. Accelerated MAP reconstructions using an accelerated factor. In: IEEE nuclear science symposium and medical imaging conference. 2010. p. 3278–81. [28] Levitan E, Herman GT. A maximum a posteriori probability expectation maximization algorithm for image reconstruction in emission tomography. IEEE Trans Med Imaging 1987;6:185–92. [29] Siddon RL. Fast calculation of the exact radiological path for a threedimensional CT array. Med Phys 1985;12:252–5. [30] Zubal IG, Harrell CR, Smith EO, Rattner Z, Gindi G, Hoffer PB. Computerized three-dimensional segmented human anatomy. Med Phys 1994;21:299–302. [31] Lin K-J, Weng Y-H, Wey S-P, Hsiao I-T, Lu C-S, Skovronsky D, et al. Whole-body biodistribution and radiation dosimetry of 18F-FP-(+)-DTBZ (18F-AV-133): a novel vesicular monoamine transporter 2 imaging agent. J Nucl Med 2010;51:1480–5. [32] Jan S, Santin G, Strul D, Staelens S, Assié K, Autret D, et al. GATE: a simulation toolkit for PET and SPECT. Phys Med Biol 2004;49:4543–61. [33] Qi J, Huesman RH. Penalized maximum-likelihood image reconstruction for lesion detection. Phys Med Biol 2006;51:4017–29. [34] Jaskowiak CJ, Bianco JA, Perlman SB, Fine JP. Influence of reconstruction iterations on 18F-FDG PET/CT standardized uptake values. J Nucl Med 2005;46:424–8. [35] Westerterp M, Pruim J, Oyen W, Hoekstra O, Paans A, Visser E, et al. Quantification of FDG PET studies using standardised uptake values in multi-centre trials: effects of image reconstruction, resolution and ROI definition parameters. Eur J Nucl Med Mol Imaging 2007;34:392–404. [36] Inoue K, Sato T, Kitamura H, Ito M, Tsunoda Y, Hirayama A, et al. Improvement of the diagnostic accuracy of lymph node metastases of colorectal cancer in 18FFDG-PET/CT by optimizing the iteration number for the image reconstruction. Ann Nucl Med 2008;22:465–73.
Please cite this article in press as: Tsai Y-J, et al. Acceleration of MAP-EM algorithm via over-relaxation. Comput Med Imaging Graph (2014), http://dx.doi.org/10.1016/j.compmedimag.2014.11.004
