ORIGINAL ARTICLE 
Annals of Nuclear Medicine Vol. 10, No.3, 299-305, 1996 
A new decision rule for parameter 8 in MAP EM (OSL) reconstruction with the Gibbs prior 
Koichi OGAWA and Keiichi HIRUMA 
Department of Electronic Informatics, College of Engineering, Hosei University 
In MAP EM (OSL) reconstruction with the Gibbs prior, the parameter delta which appears in the prior is commonly treated as a fixed value. Because the quality of reconstructed images depends on this parameter, we have to select delta very carefully, and because the statistics of an image vary locally, we should not choose a single delta value for each image. We propose a new decision rule to select an appropriate local delta. In our proposed method, delta is determined as the median of the differences between a value of the pixel of interest and those of neighboring pixels. This selection yields an appropriate prior depending on the regional statistics. The prior therefore preserves the edge property without amplifying statistical noise and it is not necessary to know the appropriate delta value to obtain high quality images. We performed computer simulations to determine the effectiveness of the proposed method. The results showed that the quality of reconstructed images obtained with the proposed method was superior to those obtained with the prior with a fixed delta. 
Key words: single photon emission CT, MAP estimation, median filter, Gibbs prior, image processing 
INTRODUCTION 
MAXIMUM a posteriori (MAP) estimation by means of an expectation maximization (EM) algorithm1 is a kind of Bayesian approach with a likelihood function and a priori information. Several prior distributions have been proposed in the literature as a priori information.2-8 The Gibbs prior9,10 is commonly used to penalize the difference between neighboring pixel values in a fixed-size neighborhood. Several methods have been proposed to perform the MAP EM reconstruction with the Gibbs prior. These include a method based on a generalized EM (GEM) algorithm4 and a method based on a one-step-late (OSL) algorithm.5-8 The MAP EM (OSL) reconstruction with an appropriate prior produces higher quality images than the maximum likelihood (ML) estimation with an EM algorithm.11,12 In the MAP EM (OSL) reconstruction with the Gibbs prior, the parameter delta in the prior is commonly treated as a fixed value.7 Because the quality of reconstructed images depends on this parameter, we must empirically select an appropriate value for parameter delta, but the statistics of an image vary locally, so it is difficult to choose an appropriate delta value. 
This paper presents a new rule for deciding the parameter delta locally in the MAP EM (OSL) reconstruction with the Gibbs prior. In the proposed method, the delta value is decided adaptively pixel-by-pixel with a median filter. The prior therefore not only preserves the edge property, but also suppresses the statistical noise. 
We performed computer simulations to confirm the validity of the proposed method. In these simulations, we compared the conventional method by using the prior with a fixed delta and the proposed method. The simulation results showed that the quality of the reconstructed images obtained with the proposed method was superior to those obtained with the conventional method. 
MAP EM (OSL) RECONSTRUCTION 
To reconstruct an image by using MAP estimation, we select an optimum image by maximizing a conditional probability of an image lambda assuming projection data Y. The posterior probability f(lambda|Y) is formulated by using Bayes' rule: 
(1) 
where f(Y|lambda) represents a likelihood function and f(lambda) represents a priori information. Since maximizing f(lambda|Y) is equivalent to maximizing In f(lambda|Y) and f(Y) is a constant for given projection data, the MAP reconstruction requires solution of the problem 
(2) 
In order to solve Eq.(2), we use the one-step-late (OSL) algorithm(5) which is an approximation to the MAP estimator. 
A. Likelihood function 
The likelihood function f(Y|lambda) reflects the statistlcal nature of the projection data Y. Since photon emission is expressed as a Poisson process, the projection data, which are a summation of the photon along a line, have a Poisson distribution. Therefore,f(Y|lambda) is given by 
(3) 
where lambda_j is the pixel value of a pixel j. Y_i is the detected number of photons for a bin i, and c_ij represents the contribution of a pixel j to a bin i.
B. Gibbs prior 
As the prior f(lambda), we assume the Gibbs prior which reflects the local statistics of an image. The Gibbs prior is defined as 
(4) 
where Z is a normalization factor which need not be calculated, B is a Gibbs distribution coefficient, and U(¥) is the energy function. The energy function is the weighted summation of the potential function V(¥), which is a function of the difference r between the values of two pixels j and l. The energy function has the following form: 
(5) 
In this equation, N_j is a neighborhood of a pixel j, w_jl is a weighting factor of the clique consisting of two pixels j and l, and delta is a scaling factor of the potential function. In our study, we used the following potential function which has been proposed by Geman and McClure:(3) 
(6) 
This function has the ability to retain the edge sharpness and reduce the statistical noise at the same time. 
METHODS 
In order to reconstruct high-quality images with the MAP EM (OSL) reconstruction, we have to select an appropriate prior area-by-area. In the MAP EM (OSL) reconstruction with the Gibbs prior, the prior has two adjustable parameters, B and delta. The Gibbs distribution coefficient B controls the contribution of all pixels neighboring an interesting pixel. The scaling factor 8 controls the contribution of a pixel neighboring the interesting pixel. In our study, we concentrated on delta because this parameter is an important factor which affects local statistics of a reconstructed image more than beta. The parameter delta is commonly treated as a fixed value in the MAP EM (OSL) reconstruction. To reconstruct high-quality images, we must empirically select an appropriate delta value. If we fix the parameter delta, we cannot realize the local statistics of an image appropriately, and, as a result, we can not preserve edges and reduce noise in the reconstructed image. 
We proposed a new decision rule for selecting the appropriate delta locally and improving the quality of the reconstructed image. In the proposed method, we used the delta value as the median of differences of two pixel values between the pixel of interest and neighboring pixels. Figure 1 shows examples of selecting delta in the different local statistics. The derivative potential function dV/dr is shown in Figure 2. In this figure, the vertical axis is the value for the derivative potential function and the horizontal axis is the difference between two pixel values (the pixel of interest and the neighboring pixels). If the pixel values are almost the same as shown in Figure 1 (a), then delta becomes small and the derivative potential function becomes the solid line in Figure 2. As a result, the flat area in the reconstructed image becomes smooth because the prior suppresses the small differences (Fig. 2, solid line r < 3). If the pixel values are distributed in Figuire 1 (b) like an edge area, then delta also becomes small. Because the prior retains the large differences (Fig. 2, solid line r > 10), the edge sharpness is preserved, but if the pixel values are distributed as shown in Figure l(c), then delta becomes large and the derivative potential function becomes the dashed line in Figure 2. As a result, we can suppress the statistical noise because the prior penalizes the large differences. In this way, we can improve the quality of the reconstructed images by using the proposed method, and it is not necessary to determine an appropriate delta to obtain a high quality image because delta is selected automatically. 
SIMULATIONS 
We performed three types of the computer simulations to determine the effectiveness of the proposed method. In these simulations, we compared the conventional method by using the prior with the fixed delta and the proposed method. 
The common simulation conditions were as follows. There were 45 views (O to 180 degrees), the initial image was formed with the mean value for the projection data, and 30 iterations were used in image reconstruction. The second order neighborhood and the weighting factor for the clique shown in Figure 3 were used in the simulations. The potential function expressed in Eq.(6) was used, and the maximum value for dV/dr was normalized to unity. In the conventional method, we used six delta values (delta=1,10,30,50,70 and 100). 
A. Simulation 1 
The computer simulations were done. by using a numerical phantom which included a flat area, an edge area and a noise area. The original 32 x 32 matrix image is shown in Figure 4. The phantom does not represent the actual activity distribution as a result of a Poisson process. Here we set typical three patterns in the phantom to basically confirm the performance of the proposed method. We examined three B values (200, 500 and 1000). 
The images reconstructed a ter rteratrons are shown in Figure 5. Figure 5(a) presents the images reconstructed by the conventional method. The delta values are, from left to right,1,30 and 100. The results of the proposed method are shown in Figure 5(c).  
The results show that the proposed method can reconstruct good images for these three areas. 
B. Simulation 2 
Three phantoms having different shapes were used for the simulations. The matrix sizes of the circular phantom and the Shepp phantom were 64 x 64, and that of the brain phantom was 128 x 128. The B value was decided empirically from a number of simulations. It was held constant at 500 in the circular phantom, 2000 in the Shepp phantom, and 2000 in the brain phantom. The brain phantom was made from an MR image. The image was quantized into four gray levels and we set arbitral values to those areas.
One percent Gaussian noise of was introduced into the projection data. The reason why we used only a tiny amount of noise was to evaluate the quality of the image easily. The original image and the images reconstructed after 30 iterations are shown in Figure 6. The original images are shown in Figure 6(a). The images reconstructed by the conventional method are shown in Figure 6(b). The three cases, from left to right, are for delta equal to 1, 30 and 100. The results obtained with the proposed method are shown in Figure 6 (c). 
The proposed method yields good results independent of the phantom shapes. 
C. Simulation 3 
The computer simulations were done with three phantoms which had different contrast but the same shape. The matrix size was 128 x 128 and the beta value was fixed at 2000 in all cases. One percent Gaussian noise was introduced into the projection data. 
The original image and the images reconstructed after 30 iterations are shown in Figure 7. Figure 7(a), (b) and (c) show the original image, the results of the conventional method, and those obtained with the proposed method. In Figure 7(b) the images reconstructed with delta = 1 , 30 and 1OO are shown at the left, middle and right. 
The proposed method produces better results than the conventional method in all cases. 
DISCUSSION 
For the images with a small delta (left in Fig. 5(a)), the edge property can be preserved in the conventional method, but the noise is not removed. This is due to the limited smoothing of only small differences. When the delta value becomes large (the middle images in Fig. 5(a)), the noise disappears, and the edge is preserved. When delta is increased further (the right images in Fig. 5(a)), the edge sharpness cannot be maintained, because the prior suppresses the large differences. When beta is large, because the reconstructed image is not affected by the prior, the flat area becomes noisy. It is therefore difficult to preserve the edge and reduce the noise at the same time with the conventional method. To sharpen up the edge without noise fluctuation in the conventional method, we have to select an appropriate delta value. In this simulation, the appropriate delta value is 30. In contrast, we can preserve the edge property and reduce the noise without respect to the beta value by using the proposed method (Fig. 5(b)) because delta is selected pixel-by-pixel according to the regional statistics of an image. 
In the proposed method, the prior at each iteration depends on the estimate of the image at that iteration, but the algorithm itself changes If the sequence of estimates changes. Consequently the uniqueness of the estimate is not preserved exactly, but simulation results show the effectiveness of the proposed method under some conditions. For the neighboring system, we used the second-order neighborhood. In the simulation we examined the first-order, the second-order, and the third-order neighborhood, and the ratio of the processing time for the selection of delta value was nearly 5 : 9 : 25. The quality of the reconstructed image, however, did not change considerably when we selected a larger order neighborhood than the second-order neighborhood. We therefore selected the second order neighborhood for the neighboring system. As for the processing time for the reconstruction, the time required was 2186 sec by the proposed method with the adaptive delta (Sun SPARCstation 10). The image size was 128 x 128, and the number of projections was 45 (30 iterations). On the other hand, the processing time for the method with the fixed delta was 1812 sec with the same conditions as in the case of the adaptive delta. For the noise treated in the simulation, we used 1% Gaussian noise for the projection data. We have examined several cases (1%,3% and more), but in the reconstructed image, it was not easy to evaluate the quality visibly without any additional image processing for the projection data except in the case of 1% noise, so we showed only the results for 1% Gaussian noise. In the clinical study, we must do some processing such as smoothing to the measured planar image before application of the proposed method. As to the convergence of the algorithm, the method may not obtaine convergence at the optimum solution. This is caused by a nonlinear process in selecting the prior, but the results shown in Figure 8 imply the stability of the algorithm under the conditions used in the paper. 
As shown in Figure 6(b), a small delta amplified noise and a large delta reduced the edge sharpness of reconstructed images. We must therefore select an appropriate delta carefully for each phantom in the conventional method because the reconstructed image is sensitive to the delta value. Nevertheless, the proposed method yields highquality images without information on the appropriate delta value for an shape of phantom or for an distribution. 
Figure 7 shows that the proposed method produces better results than the conventional method, as in simulations 1 and 2. To illustrate this further, we calculated the mean square error between the original image and the reconstructed images. Figure 8 shows the mean square error for each phantom in Figure 7. The pixel values are 140, 150 and 160 in (a); 100, 150 and 200 in (b); and 50, 150 and 200 in (c). The vertical axis in Figure 8 represents the mean square error, and the horizontal axis, the number of iterations. As shown in Figures 8(a) and (b), both the mean square error of the conventional method and the proposed method are stable for the iterations. The mean square error of the images when using the proposed method is smaller than that when using the conventional method. In Figure 8(c), the mean square error of the conventional method begins to increase after a certain number of iterations. The instability in the mean square error is caused by the amplification noise. In order to obtain good results with the conventional method, we need to know the appropriate number of iterations. In contrast, because the mean square error of the proposed method is stable, the proposed method yields high-quality images without requiring knowledge of the appropriate number of iterations. 
CONCLUSIONS 
In this paper, we proposed a new rule for deciding the parameter delta in MAP EM (OSL) reconstruction. The new rule uses the Gibbs prior to avoid selecting delta empirically and improve the quality of reconstructed images. In the simulations, the proposed method produced better results than the conventional method with the prior with the fixed delta because it selects delta pixel-by-pixel depending on the local statistics of an image. Moreover, it is not necessary to decide an appropriate delta value to reconstruct high-quality images. Further work will be needed to estimate the appropriate value for the Gibbs distribution coefficient B. 
REFERENCES 
l . Dempster AP, Laird NM, Rubin DB. Maximum likelihood from incomplete data via the EM algorithm (with discussion. J Royal Statist Soc B 39: 1-38, 1977. 
2. Lange K, Bahn M, Little R. A theoretical study of some maximum likelihood algorithms for emission and transmission tomography. IEEE Trans Med Imag MI-6:106-114, 1987. 
3. Geman S, McClure DE. Bayesian image analysis: An application to single photon emission tomography, Proceedings of the Statistical Computing Section, Washington, DC., Amer Statist Assoc, pp. 12-18, 1985. 
4. Hebert TJ, Leahy R. A generalized EM algorithm for 3-D Bayesian reconstruction from Poisson data using Gibbs priors. IEEE Trans Med Imag MI-8: 194-202, 1989. 
5. Green PJ. On the use of the EM algorithm for penalized likelihood estimation. J Royal Statist Soc B 52: 443-452, 1990. 
6. Lange K. Convergence of EM image reconstruction algorithms with Gibbs smoothing. IEEE Trans Med Imag MI-9: 439-446, 1990. 
7. Lalush DS. Tsui BMW. Simulation evaluation of Gibbs prior distributions for use in maximum a posteriori SPECT reconstructions. IEEE Trans Med Imag MI-11 : 267-275, 1992. 
8. Lalush DS, Tsui BMW. A generalized Gibbs prior for maximum a posteriori reconstruction in SPECT. Phy Med Biol 38: 729-741 , 1992. 
9. Geman S, Geman D. Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images. IEEE Trans Pattern Anal Mach Intell PAMI-6: 721-741 , 1984. 
10. Besag J. On the statistical analysis of dirty pictures (with discussion). J Royal Statist Soc B 48: 259-302, 1986. 
11. Shepp LA, Vardi Y. Maximum likelihood reconstruction for emission tomography. IEEE Trans Med Imag MI-1 : 113-122, 1982. 
12. Lange K, Carson R. EM reconstruction algorithms for emission and transmission tomography. J Comput Assist Tomog 8: 306-316, 1984.