Nonlocal Prior Bayesian Tomographic ... - Semantic Scholar

J Math Imaging Vis (2008) 30: 133–146 DOI 10.1007/s10851-007-0042-5

Nonlocal Prior Bayesian Tomographic Reconstruction Yang Chen · Jianhua Ma · Qianjin Feng · Limin Luo · Pengcheng Shi · Wufan Chen

Published online: 29 November 2007 © Springer Science+Business Media, LLC 2007

Abstract Bayesian approaches, or maximum a posteriori (MAP) methods, are effective in providing solutions to illposed problems in image reconstruction. Based on Bayesian theory, prior information of the target image is imposed on image reconstruction to suppress noise. Conventionally, the information in most of prior models comes from weighted differences between pixel intensities within a small local neighborhood. In this paper, we propose a novel nonlocal prior such that differences are computed over a broader neighborhoods of each pixel with weights depending on its similarity with respect to the other pixels. In such a way connectivity and continuity of the image is exploited. A twostep reconstruction algorithm using the nonlocal prior is developed. The proposed nonlocal prior Bayesian reconstruction algorithm has been applied to emission tomographic reconstructions using both computer simulated data and patient SPECT data. Compared to several existing reconstruction methods, our approach shows better performance in both lowering the noise and preserving the edges. Keywords Bayesian reconstruction · Positron emission tomography (PET) · Single photon emission computed tomography (SPECT) · Local prior · Nonlocal prior

Y. Chen · L. Luo Laboratory of Image Science and Technology, Southeast University, Nanjing 210096, China Y. Chen · J. Ma · Q. Feng · P. Shi · W. Chen () School of Biomedical Engineering, Southern Medical University, Guangzhou 510515, China e-mail: [email protected]

1 Introduction Because of low counting rates and noise, the reconstruction of an unknown image f from measurement data g, such as the reconstruction of Positron Emission Tomography (PET) images and Single Positron Emission Computerized Tomography (SPECT) images, is an ill-posed problem [1, 2]. The maximum-likelihood expectation-maximization (ML-EM) approach, better expressing system models of physical effects and modeling the statistical character of the measurement data, outperforms conventional analytic methods such as the filtered backprojection (FBP) method [3]. However, the ML-EM approach is still notorious for its slow convergence, and the reconstructed images often start to deteriorate, producing “checkerboard effect”, as the iteration proceeds [4]. Overcoming such ill-posedness for image reconstruction has always been a problem. Bayesian methods, or maximum a posteriori (MAP) methods, have long been accepted as effective solutions to this problem [4–6]. Based on Bayesian and MAP theory, the prior information can be transformed into constraints to regularize the solution to the original illposed reconstruction problem. The prior information measures the extent to which the constraints are violated by the image. Therefore, regularization from such prior information can be imposed on image reconstruction to suppress noise. One can build the following posterior probability P (f |g) for image reconstruction: P (f |g) ∝ P (g|f )P (f ),

(1) P (f ) = Z −1 × e−βU (f ) = Z −1 × exp −β Uj (f ) (2) j

where g and f denote the measurement data and the image to be reconstructed, respectively. P (g|f ) and P (f ) are

134

J Math Imaging Vis (2008) 30: 133–146

the likelihood distribution and the prior distribution, respectively. Z is a normalizing constant called the partition function. Uj (f ) is the notation for the value of the energy function U (f ) evaluated at pixel index j , and β is the global hyperparameter that controls the degree of the prior’s influence on the image f . Then, we can build the posterior energy function as ψβ (f/g) = log P (f/g) = L(g, f ) − βU (f )

(3)

here L(g, f ) represents the likelihood energy function. Then, the reconstructed image f can be obtained through maximization of function ψβ (f ) by an iterative procedure. The posterior objective energy function (3) does not always produce satisfactory results because many images are not globally smooth, and edges combined with noise often make the values of some region vary abruptly. The quadratic membrane (QM) smoothing prior, which smoothes both noise and edge details equally, tends to produce an unfavorable oversmoothing effect. On the other hand, many edgepreserving Bayesian methods, including the line-site model [6, 7] and DA (Discontinuity Adaptivity) MRF model [6], can produce sharper edges using a nonquadratic prior [6– 12]. These edge-preserving Bayesian methods with nonquadratic priors might work well when the noise level of detected data is low and the image to be reconstructed is composed of homogeneous regions with sharp boundaries, otherwise they tend to produce blocky piecewise smooth regions or so-called staircase effect. In clinical situation, this staircase effect in medical imaging might be misinterpreted. In addition, the edge-preserving Bayesian methods are usually computationally costly because of often needed annealing procedures and estimations of some built-in parameters of priors [9–11]. The quadratic prior can smooth the reconstructed image by averaging pixels within a local neighborhood. The edgepreserving nonquadratic priors rely on pixel density difference information within a local neighborhood to determine the degrees of regularization on every pixel in the image [6, 13]. However, none of them addresses global connectivity and continuity information in the image. Both the oversmoothing effect caused by the quadratic priors and the staircase effect caused by the nonquadratic priors are the outcomes of their limited power of distinguishing edge from noise. Such limited power can be ascribed to the fact that the simple weighting strategy for pixel intensity differences within a small fixed local neighborhood can only provide limited prior information. In this paper we call above conventional priors, which depend on the simple weighting strategy within a local neighborhood, local priors. Our aim is to find a new prior that is capable of exploiting more information in the images. There exist several methods which claim to capture global properties of the

images, including the methods using region-based Bayesian prior [14, 15] and boundary-based Bayesian prior or regularization that incorporates global image information by levelsets [13]. However, the region-based methods are restricted by the complicated region identification and anatomical prior information. On the other hand, the boundary-based method relies heavily on the level-set operations whose effect in different images is unpredictable and parametersensitive [13]. Recently, Buades and his colleagues put forward a novel nonlocal algorithm for image denoising [16, 17]. Inspired by their nonlocal idea and based on their work, we build a nonlocal prior model for Bayesian image reconstruction. Exploiting not only intensity-differences information between individual pixels but also global connectivity and continuity information within a relatively larger neighborhood in the image f , the nonlocal prior can greatly improve reconstruction. In Sect. 2, the theory of the conventional local prior model and the theory for the proposed nonlocal prior model are presented. Reconstruction algorithm for the Bayesian reconstruction using the nonlocal prior is developed in Sect. 3. In Sect. 4, we reconstruct images using the proposed algorithm using both simulated data and real SPECT data. Qualitative and quantitative comparisons and analyses are also included. Conclusions and plans for future work are given in Sect. 5.

2 Nonlocal Prior Model Before introducing the nonlocal prior model, we review the conventional local prior model. 2.1 Conventional Local Prior Model According to Bayesian theory, when the image f meets the prior assumptions, the prior energy function U (f ) in (2) attains its minimum and the corresponding prior distribution (2) attains its maximum. Conventionally, the value of Uj (f ) is commonly computed through a weighted sum of potential functions v of the differences between pixels in the neighborhood Nj : Uj (f ) = β

wbj v(fb − fj ).

(4)

b∈Nj

Generally, different choices of the potential function v lead to different priors. The prior becomes the conventional space-invariant QM prior when the potential function takes the form of v(t) = t 2 . We can also choose edge-preserving nonquadratic priors by adopting a nonquadratic potential


135

function for v, such as the well-known Huber potential function: |t| ≤ δ, t 2 /2, (5) v(t) = 2 δ|t| − δ /2, |t| > δ where δ is the threshold parameter [4–10, 13]. Weight wbj is a positive value that denotes the interaction degree between pixel b and pixel j . In conventional local prior model, it is usually considered to be inversely proportional to the distance between pixel b and pixel j . On a square grid of image f , in which the two dimensional positions of the pixel b and pixel j (b = j ) are respectively (bx , by ) and (jx , jy ), wbj are usually calculated simply by wbj = 1/dbj = 1/ (bx − jx )2 + (by − jy )2 .

(6)

A typical normalized eight-neighborhood weighting map and a typical four-neighborhood weighting map for local priors are as follows:

(1)

(2)

√ √ 1/ 2 1 1/ 2 √ (1/(4 × 1 + 4 × 1/ 2)) × 1√ 0 1√ , 1/ 2 1 1/ 2 1 (1/(4 × 1)) × 1 0 1 . 1

These kinds of simple weight-determining methods have been used widely in Bayesian image reconstructions and restorations [9]. From above we can see that the local priors can only provide fixed local prior information for image reconstruction. The local QM prior is prone to producing oversmoothing effect which smoothes out both noise and the important details [6, 12]. The local edge-preserving nonquadratic priors, while able to preserve some edge information by choosing a nonquadratic potential function that increases less rapidly when the differences between adjacent pixels become bigger, tend to produce some unfavorable artifacts [13]. 2.2 Proposed Nonlocal Prior Model Recently, Buades and his colleagues applied a novel nonlocal algorithm in image denoising [16, 17]. Inspired by their nonlocal idea, we devise a novel nonlocal prior for Bayesian image reconstruction. In the development of the nonlocal prior, a large neighborhood N is set to incorporate geometrical configuration information in the image. And rather than just using above mentioned simple reciprocal of the distance between two individual pixels, the weight wbj for the nonlocal prior is set to reflect the degree of similarity or connectivity between pixel b and pixel j in neighborhood N

base on the idea in [16, 17]. So, we define wbj as a decreasing function of the Euclidean distance of the two square comparing neighborhoods nb and nj centered at pixel b and pixel j , respectively. The two neighborhoods nb and nj have the same sizes. The building of the proposed nonlocal prior depends on image f in some way and can be formulized as follows: UN L (f ) = Uj (f ) = wbj (fb − fj )2 , (7) j

j b∈Nj

f (nb ) − f (nj )2E dbj , wbj = exp − h2

(8)

f (nb ) = {fl : l ∈ nb },

(9)

f (nj ) = {fl : l ∈ nj },

(10)

(fl∈nb − fl∈nj )2 . f (nb ) − f (nj )2E =

(11)

l

Here, dbj can be computed using (6). f (nb ) and f (nj ) denote the two pixel intensity vectors in the two comparing neighborhoods nb and nj , respectively. Then wbj , which represents the similarity between the two neighborhoods nb and nj , can be computed via (8–11). Parameter h in (8) controls the decay of the exponential function of the weights for the nonlocal prior. Through the computation of the similarities of the comparing neighborhoods n centered on the pairs of pixels in the large neighborhood N , the spatial-varying weights of each pair of pixels in N for the proposed nonlocal prior model are distributed across the more similar configurations. The proposed nonlocal prior is able to systematically use global self-prediction configuration information in the image. To make above statement more lucid and straightforward, we illustrate several computed prior weight distribution maps for tagged pixels in Fig. 1. In all the illustrations, the sizes of neighborhoods N and n are set to be 21×21 and 7×7, respectively. In Fig. 1(a), we can see that the prior weight distributions address the direction of the straight edge line in neighborhood N when the pixel is in a straight edge. In Fig. 1(b), when the pixel belongs to a curved edge, the corresponding prior weight distributions favor pixels along the same curved line in neighborhood N . In Fig. 1(c), the prior weight distributions show no configuration information in neighborhood N when the pixel is located in region with no configuration. Figure 1(d) illustrates two prior weight distributions for two pixels in edge region and background region, respectively. We can also see that the prior weights for the pixel in edge region have large values along the edge region in neighborhood N , while the prior weights for the pixel in background region do not show any such pixel configuration representation.

136


Fig. 1 Illustrations of the prior weight distributions for the pixels (the white points pointed by arrows) using the proposed nonlocal prior model. a Prior weight distributions (the right image) in neighborhood N for the pixel in a straight edge. b Prior weight distributions (the right image) in neighborhood N for the pixel in a curved edge. c Prior

weight distributions (the right image) in neighborhood N when the pixel is not located in any configuration region. d Prior weight distributions (the upper-left and the lower-right images) in neighborhood N for the two pixels in edge region and background region, respectively

We must manifest that the nonlocal character of the method does not mean or necessitate unbounded neighborhoods. Its motive or essence is to incorporate more connectivity and continuity information into the prior for reconstructions. For computation complexity consideration, we generally do not extend the bounds of neighborhood N and n over the whole images, and we set the sizes to be appropriate values in study. In experiments, we find that the 11×11 neighborhood N and the 7 × 7 comparing neighborhood n have shown to be large enough to be robust to noise.

the measurement vector g with the emission isotope distribution imagef . Therefore, following formulas can be deduced:

3 Two-Step Joint Reconstruction Algorithm 3.1 Statistical Model It is known that emission tomography is usually an ill-posed inverse problem due to low count rates, noise contamination and other physical effects. In statistical reconstruction, PET or SPECT measurements are based on a counting process, and it is reasonable to assume that the measurement vectors are independently distributed Poisson random variables. Let gi be the measurement detected by the detector pair i in the case of PET or detector i in the case of SPECT. It can be well modeled as Poisson random with expectation g¯ i as a function of the underlying isotope distribution image vector [3]. In emission tomography, PL (g/f ) (i.e., the likelihood function) is the conditional probability of obtaining

gi ∝ Poisson{g¯ i (f )}, g¯ i (f ) =

J

i = 1 . . . I,

ci aij fj + ri ,

(12) (13)

j

PL (g/f ) =

I g¯ i (f )gi i=1

gi !

exp(−g¯ i (f ))

(14)

where N is the number of detector pairs, M is the number of image pixels, aij is the geometric probability that an emission photon from image pixel j is detected by the detector pair i, ci incorporates the effect of scan time, detector efficiencies, attenuation factors and possibly deadtime correction factors for detector pair i, and ri represents the total detected random and scattered event counts for detector pair i, respectively. Then, the corresponding log-likelihood function for estimating f from g is [3]: L(f ) = log PL (g/f ) =

I (gi log g¯ i (f ) − g¯ i (f ) − log gi !).

(15)

i

3.2 Joint Reconstruction Algorithm Choosing the nonlocal prior, the objective image f can be obtained through an iterative maximization of the posterior


137

energy function ψβ (f/g): f = arg max ψβ (f/g) f

= arg max(L(g, f ) − βUN L (f )).

(16)

gence proofs [4, 24, 25]. Although no guarantee of convergence, the proposed joint updating reconstruction algorithm is, however, effective in practice and able to converge to at least a local maximum as the joint algorithm proposed in [13].

f

However, from (7) and (8), we find that every weighting factor wbj is determined by the imagef , which makes it difficult to obtain partial derivatives. We apply following alternating two-step updating reconstruction algorithm in our study. After choosing an initial image estimate f 0 , we update weight w and image f alternatively by following two-step strategy for every iteration. (1) Weight update For every pixel pair (fb , fj ) in image f with fˆ being the current estimate of f , compute wbj using (7–11) and current image estimate fˆ. (2) Image update It has been shown in [3] that ∂ 2 L(f )/∂fb ∂fj with L(f ) given by (14) is nonpositive definite. From (8–11), we know for every pixel pair (fb , fj ) 0 < wbj (fˆ) fˆ(nb ) − fˆ(nj )2 1 dbj ≤ = exp − . dbj h2

(17)

Considering the convexity for (fb − fj )2 in (7) and the positivity for hyperparameter β, we conclude that, given the computed value for wbj (fˆ), the second derivatives of ψβ (f/g) (16) is definitely concave in each iteration. Thus based on the theory of locally linearization [18], we can monotonically maximize the posterior energy function ψβ (f/g) using iterative algorithms such as the paraboloidal surrogate coordinate ascent (PSCA) iterative algorithm [19, 20] or the conjugate gradient iterative algorithm [21–23]. We chose the PSCA algorithm in our study. This joint updating approach, which can be considered as a version of the widely-used OSL (One Step Late) iterative algorithm, suffers from the lack of strict converFig. 2 The synthetic phantoms used in computer simulations: a Phantom 1, b Phantom 2, and c Phantom 3

4 Experimentation Both simulated emission data and real emission data are used in our study to evaluate the proposed nonlocal prior Bayesian reconstruction algorithm. 4.1 Emission Reconstruction Using Simulated Data In this experiment, three synthetic simulated phantoms are used for PET reconstruction. Figure 1 shows the three computer simulated phantoms. Figure 2(a) is a Zubal sectional phantom with pixel values from 0 to 8, and the total counts amount to 4 × 105 ; Fig. 2(b) is a phantom composed by a constant background region, a square hot region and a blobby hot region containing another small square cold region. This phantom has pixel values ranging from 0 to 5, and the total counts amount to 4 × 105 ; Fig. 2(c) is a SheppLogan head phantom with pixel values from 0 to 8, and the total counts amount to 3 × 105 . For simplicity, we refer to above three phantoms as phantom 1, phantom 2 and phantom 3, respectively. The simulated sonogram data for above three phantom data are all Poisson distributed and the percentages of simulated delayed coincidences ri factors (scatter effects are ignored) in all counts are all set to be 10%. All ci ψ factors are generated using pseudo-random log-normal variates with a standard deviation of 0.3 to account for detector efficiency variations. Generated by the ASPIRE software system [26], the three transition probability matrices used in the reconstructions of above phantoms all correspond to parallel stripintegral geometry with 128 radial samples and 128 angular samples distributed uniformly over 180 degrees. Images are stored in 128 × 128 matrices. A PC with Intel P4 CPU and 512 Mb RAM is used as the workstation for reconstruction. Same reconstruction software environments are set for the reconstructions of above three phantoms.

138


Table 1 CPU times needed for the simulated phantom data reconstructions using different reconstruction methods

CPU time

FBP recon.

QM prior recon.

Huber prior recon.

Nolocal prior recon.

10.81 s

172.55 s

250.23 s

612.04 s

The FBP reconstructions and Bayesian reconstructions using the QM prior, Huber prior, and the proposed nonlocal prior are performed. For the FBP method, we choose the ramp filter with the cutoff frequency equal to the Nyquist frequency. Different parameters are used for the reconstructions for different phantoms. For the study using the same phantom, same hyperparameter β can be used in reconstructions using the QM prior and Huber prior. Hyperparameter β and threshold parameter δ for reconstructions using the QM prior and Huber prior are chosen by hand to give the reconstructions with the maximum signal-to-noise ratio (SNR), which is computed by: J SNR = 10 log10

2 ¯ j (fphantom−j − fphantom ) . J 2 j (fj − fphantom−j )

(18)

In (18), the numerator and denominator stand for the signal power and the noise power, respectively. f, f¯, fphantom and f¯phantom denote the reconstructed image, the mean of all pixels in f , the corresponding true phantom images (Fig. 2), and the mean of all pixels in the phantom images, respectively. For all reconstructions using the QM prior and Huber prior, the eight-neighborhood mentioned in Sect. 2.1 is used, and the value of β is set to 1.5, 2 and 1.75 for the studies using the three phantom, respectively. The value of δ for Huber prior is fixed to 0.2 in all reconstructions. As to the Bayesian reconstructions using the proposed nonlocal prior, the values of parameters are also chosen empirically to give the most stable reconstructed images with the maximum SNR. The value of parameter β is set to 0.5, 0.6 and 0.4 for the studies using the three phantoms, respectively. The value of parameter h of the nonlocal prior in (8) is set to be 0.6, 1.1 and 0.8 for the studies using the three phantom, respectively. For all Bayesian reconstructions using the nonlocal prior, the neighborhood N in (7) is set to an 11×11 neighborhood and the two comparing neighborhoods nb and nj in (8) are both set to a 7×7 neighborhood. It is found in practice that the improvement of reconstruction is not obvious when the sizes of neighborhood N and comparing neighborhood n are larger than these size settings. The reconstructed images from the FBP method are used as the initial images in iteration processes, and the 150th iterated images in reconstruction are chosen as the results. After 150 iterations, the images do not change visually. Above joint reconstruction algorithm is used in experiments.

Figure 3(1), (2), and (3) show the reconstructed images for the three phantoms, respectively. Each image (a) in (1), (2) and (3) is the reconstructed image using the FBP algorithm. Images (b), (c), (d) in every row are the Bayesian reconstructions using the QM prior, Huber prior, and the proposed nonlocal prior, respectively. Considering the same sizes of phantom data and system matrices, same computational cost is needed for the reconstructions of different phantom data when same type of prior is used. So we only give one CPU time cost for different phantom data when using same reconstruction method. The corresponding CPU times for all reconstructions in MATLAB 7.0 software environment are displayed in Table 1. We can see that the computational costs for Bayesian reconstructions using the proposed nonlocal prior are comparably higher. 4.1.1 Qualitative Comparisons From Fig. 3, we can see that the FBP reconstructions suffer severely from the noise and streak artifacts, and the reconstructions using the proposed nonlocal prior yield more appealing images than the other methods. Obviously, reconstructions using the proposed nonlocal prior exhibit excellent performance in both suppressing the noise and preserving the edges, and are free from the unfavorable oversmoothing effect which is usually seen with the QM prior and the staircase effect with nonquadratic Huber prior. Figure 4 depicts the horizontal profiles of the Bayesian reconstructions in Fig. 3. From the profiles, we can see that the resulting images from the reconstructions using the nonlocal prior have the profiles that agree better with the profiles of the phantoms at the upper left corner in Fig. 4(1), (2), and (3), respectively. 4.1.2 Quantitative Comparisons Quantitative comparisons [27] are also made for further analyses of reconstructions using the proposed nonlocal prior. Table 2 displays SNR comparisons and indicates that the reconstructed images using the proposed nonlocal prior have higher SNRs than other reconstruction methods. We choose Phantom 2 to test the bias, average standard deviation (astd) and contrast trade-offs for reconstructions using the three priors (QM prior, Huber prior and the nonlocal prior). As is shown in Fig. 5, a region of interest (ROI) and a background region, which are surrounded by solid


139

Fig. 3 Emission FBP reconstruction and Bayesian reconstructions using different priors in the three phantom studies: a FBP reconstruction, b QM prior reconstruction, c Nonquadratic Huber prior reconstruction, and d the proposed nonlocal prior reconstruction

white line, are chosen for the study. Fifty pseudo-random Poisson-distributed sinogram sets are simulated for each β. Let W denote the chosen ROI or background region, f e denote the 150th iterated reconstruction for the eth realization data set, f¯ denote the sample mean of all the realization, fphantom denote the chosen phantom, C = 50 be the number of pseudo-random sinogram realizations, NW be the number of pixels in the ROI or the background region, and MeanROI and MeanBackground denote the mean of all the pixels in the ROI and the background region, respectively. Then:

contrast =

|MeanROI − MeanBackground | , MeanBackground

(19)

|f¯j − f¯phantom−j | , ¯ j ∈W |fphantom−j |

C

1 1 = (f¯j − fje )2 . NW C −1

bias(f )W =

astd(f )W

j ∈W

j ∈W

(20)

(21)

e=1

Comparisons are made for different h for the nonlocal prior and different δ for Huber prior. Hyperparameter β is used to regulate the values of biases, astd and contrast. Figure 6(a) shows the legend for all the trade-off plots in Fig. 6. Figure 6(b) and (c) are the plots of the astd versus bias for the ROI and the background region, respectively. We can see that reconstructions using the pro-

140


1) Horizontal profiles of the images from Bayesian reconstruction in phantom study 1

2) Horizontal profiles of the images from Bayesian reconstruction in phantom study 2 Fig. 4 Horizontal profiles of the reconstructed images from Bayesian reconstructions in Fig. 3


141

3) Horizontal profiles of the images from Bayesian reconstruction in phantom study 3 Fig. 4 (Continued) Table 2 SNRs for the reconstructed images Fig. 2

SNR

FBP recon.

QM prior recon.

Huber prior recon.


Phantom 1

11.34

12.44

13.50

14.10

Phantom 2

11.67

13.07

14.12

15.23

Phantom 3

10.72

12.81

13.21

14.85

Fig. 5 The chosen ROI and Background region (the areas surrounded by solid white line) in Phantom 2: a ROI; b Background region

posed nonlocal prior can obtain both lower astd and bias, and show better performance in terms of astd versus bias

trade-off than reconstructions using the QM prior and Huber prior. Figure 6(d) depicts the plots of the astd in the ROI versus the astd in the background region. In Fig. 6(d), we can see that the curves for reconstructions using the nonlocal prior lie lower-left to the curves for reconstructions using the QM prior and Huber prior. So, we conclude that reconstructions using the proposed nonlocal prior outperform reconstructions using the QM prior and Huber prior in terms of the trade-off of the astd in the ROI versus the astd in the background region. For the trade-off of the contrast between the ROI and the background region versus the astd of the background region, it can be clearly seen in Fig. 6(e) that the trade-off curves for reconstructions using nonlocal prior lie above the curves for

142

Fig. 6 Plots of the bias, astd and contrast trade-offs for reconstructions using the QM prior, Huber prior and the proposed nonlocal prior (the values of hyperparameter β for the nonlocal prior, the QM prior and Huber prior range from 0.015 to 2, 0.025 to 3 and 0.015 to 1.8, respectively). a The legend for the plots (b), (c), (d), and (e); b Plots of


the astd in the ROI versus the bias in the ROI; c Plots of the astd in the background region versus the bias in the background region; d Plots of the astd in the ROI region versus the astd in the background region; e Plots of the contrast between the ROI and the background region versus the astd in the background region


reconstructions using the QM prior and Huber prior, which indicates that reconstructions using the proposed nonlocal prior are also preferable in terms of the contrast trade-off and able to reconstruct images with better contrast and distinguishability. 4.2 Reconstruction Using Real Emission Scan Data In this study, real patient SPECT data obtained from cardiac perfusion imaging using Tc-99m sestamibi are used. For each slice, 128 projections were taken at 120 uniformly spaced angles between 0 and 2π . The 9th 2D slice data of the total 39 slices is used. The image consists of 128 × 128 pixels with 0.356 cm pixel resolution. The total counts amount to 1.4876 × 105 . All ci factors and ri factors are assumed to be ones and zeros, respectively. For all Bayesian reconstructions, the reconstructed images from the same FBP method as described above are used as the initial image in the iterative algorithms. The QM prior, Huber prior and the nonlocal prior for the patient data are the same as those used in the simulated data study. Iteration numbers for all Bayesian reconstructions are set to be 50. The proposed joint reconstruction algorithm is also used here. Figure 7 shows the reconstructed results from different methods with different parameter settings. For all Bayesian

1) FBP reconstruction with ramp filter with cutoff frequency equal to the Nyquist frequency Fig. 7 FBP reconstruction (the first row) and Bayesian reconstructions using different priors under different parameter levels (the second, third, fourth, fifth and sixth rows) for real emission data: The second row are the Bayesian reconstructions using the QM prior for β = 0, 25, 50, 70 (from left to right); The third row are the Bayesian reconstructions using Huber prior for δ = 0.02 and β = 5, 10, 25, 35 (from left to right); The fourth row are the Bayesian reconstructions using Huber prior for β = 50 and δ = 0.005, 0.01, 0.02, 0.03 (from left to right); The fifth row are the Bayesian reconstructions using the nonlocal prior for h = 0.02 and β = 5, 10, 25, 35 (from left to right); The sixth row are the Bayesian reconstructions using the nonlocal prior for β = 25 and h = 0.005, 0.01, 0.02, 0.03 (from left to right)

143

reconstructions, we manually select parameters β, δ and h to cover a range of parameter settings. Obviously, when suitable parameter values are set, the reconstructions using the proposed nonlocal prior outperform FBP reconstruction and reconstructions using QM prior and Huber prior in terms of noise suppression and edge preserving, and are able to reconstruct more appealing emission images. The corresponding CPU times for the different Bayesian reconstructions in MATLAB 7.0 software environment are displayed in Table 3. We can also see that the computation cost for reconstructions using the proposed nonlocal prior is comparably higher.

5 Conclusion Stemming from Bayesian theory and the approach in [16, 17], the proposed nonlocal prior is theoretically reasonable and straight forward. Through choosing a large neighborhood N and a new weighting strategy that incorporates the connectivity between pixels in N , the proposed nonlocal prior is able to provide more image global information for the original ill-posed reconstruction. Tomographic reconstruction using the nonlocal prior can be feasibly and effectively implemented using the proposed two step joint reconstruction algorithm in Sect. 3.2. And, from the analyses and experiments presented in this paper, our proposed nonlocal prior demonstrates a more effective and robust regularization for PET and SPECT reconstruction than the conventional local QM prior and Huber prior. Nevertheless, as can be seen from Table 1 and Table 3, computation over a large neighborhood and the incorporation of the new empirically adjusted parameter h increases the computational complexity for the nonlocal prior reconstructions. However, considering the fast development of computer technology and the urgent clinical needs for highquality PET and SPECT images, the increased computational cost should not be overblamed. Ordered subset (OS) [28] and space-alternating generalized EM (SAGE) [29] accelerating methods can be combined with the reconstruction algorithm to accelerate the convergence rates and reduce the computational cost. Future work includes analyzing the convergence property of the reconstruction algorithm, exploring effective ways to reduce the computational cost and studying automatic parameter estimation methods. We will also make more experiments to test the nonlocal approach’s applications in transmission tomography, and image deconvolution in microscopy and astronomy.

144

Fig. 7 (Continued)



145

6) Bayesian reconstructions using the proposed nonlocal prior for β = 25 and h = 0.005, 0.01, 0.02, 0.03, respectively (from left to right) Fig. 7 (Continued) Table 3 CPU times needed for real emission data reconstructions using different reconstruction methods

CPU time

FBP recon.

QM prior recon.

Huber prior recon.


9.59 s

74.78 s

90.65 s

222.63 s

Acknowledgements This research was supported by National Basic Research Program of China under grant, No. 2003CB716102. Authors make a grateful acknowledgement for revision work of Professor Larry Zeng’s revision work.

References 1. Bertero, M., De Mol, C., Pike, E.R.: Linear inverse problems with discrete data, I: general formulation and singular system analysis. Inverse Prob. 1(4), 301–330 (1985) 2. Bertero, M., Poggio, T., Torre, V.: Ill posed problems in early vision. Proc. IEEE 76, 869–889 (1988) 3. Shepp, L.A., Vardi, Y.: Maximum likehood reconstruction for emission tomography. IEEE Trans. Med. Imaging 1, 113–121 (1982) 4. Lange, K.: Convergence of EM image reconstruction algorithms with Gibbs smoothness. IEEE Trans. Med. Imaging 9, 439–446 (1990) 5. Geman, S., Geman, D.: Stochastic relaxation, Gibbs distribution, and the Bayesian restoration of images. IEEE Trans. Pattern Anal. Mach. Intell. 6, 721–741 (1984) 6. Li, S.Z.: Markov Random Field Modeling in Image Analysis. Springer, Tokyo (2001) 7. Black, M.J.: Rangarajan A., Unification of line process, outlier rejection, and robust statistics with application in early vision. Int. J. Comput. Vis. 19, 57–91 (1996) 8. Charbonnier, P., Blanc-Féraud, L., Aubert, G., Barlaud, M.: Deterministic edge preserving regularization in computer imaging. IEEE Trans. Image Process. 6(2), 298–311 (1997) 9. Johnson, V., Wong, W.H., Hu, X., Chen, C.T.: Image restoration using Gibbs prior: boundary modeling, treatment of blurring, and selection of hyperparameter. IEEE Trans. Pattern Anal. Mach. Intell. 13, 413–425 (1991) 10. Bouman, C.A., Sauer, K.: A generalized Gaussian image model for edge-preserving MAP estimation. IEEE Trans. Image Process. 2, 296–310 (1993) 11. Gindi, G., Rangarajan, A., Lee, M., Hong, P.J., Zubal, G.: Bayesian reconstruction for emission tomography via deterministic annealing. In: Barrett, H., Gmitro, A. (eds.) Information Processing in Medical Imaging, pp. 322–338. Springer, Berlin (1993)

12. Lee, S.J., Rangarajan, A., Gindi, G.: Bayesian image reconstruction in SPECT using higher order mechanical models as priors. IEEE Trans. Med. Imaging 14(4), 669–680 (1995) 13. Yu, D.F., Fessler, J.A.: Edge-preserving tomographic reconstruction with nonlocal regularization. IEEE Trans. Med. Imaging 21(2), 159–173 (2002) 14. Johnson, V.E.: A model for segmentation and analysis of noisy images. J. Am. Stat. Assoc. 89(425), 230–241 (1994) 15. Bowsher, J.E., Johnson, V.E., Turkington, T.G., Jaszczak, R.J., Floyd, C.E., Coleman, R.E.: Bayesian reconstruction and use of anatomical a priori information for emission tomography. IEEE Trans. Med. Imaging 15(5), 673–686 (1996) 16. Buades, A., Coll, B., Morel, J.M.: A nonlocal algorithm for image denoising. In: Proc. IEEE Int. Conf. Computer Vision Pattern Recognition, vol. 2, pp. 60–65 (2005) 17. Buades, A., Coll, B., Morel, J.M.: A review of image denoising algorithms, with a new one. Multiscale Model. Simul. 4(2), 490– 530 (2005) 18. Chen, W., Chen, M., Zhou, J.: Adaptively regularized constrained total least-squares Image Restoration. IEEE Trans. Image Process. 9(4), 588–596 (2000) 19. Fessler, J.A., Erdo˘gan, H.: A paraboloidal surrogates algorithm for convergent penalized-likelihood emission reconstruction. In: Proc. IEEE Nucl. Sci. Symp. Med. Imag. Conf., vol. 2, pp. 1132–1135 (1998) 20. Erdo˘gan, H., Fessler, J.A.: Monotonic algorithms for transmission tomography. IEEE Trans. Med. Imaging 18(9), 801–814 (1999) 21. Lalush, D.S., Tsui, B.M.W.: A fast and stable maximum a posteriori conjugate gradient reconstruction algorithm. Med. Phys. 22(8), 1273–1284 (1995) 22. Fessler, J.A., Booth, S.D.: Conjugate-gradient preconditioning methods for shift-variant PET image reconstruction. IEEE Trans. Image Process. 8(5), 688–699 (1999) 23. Booth, S.D., Fessler, J.A.: Combined diagonal/Fourier preconditioning methods for image reconstruction in emission tomography. In: Proc. IEEE Int. Conf. on Image Processing, vol. 2, pp. 441–444 (1995) 24. Alenius, S., Ruotsalainen, U.: Generalization of median root prior reconstruction. IEEE Trans. Med. Imaging 21(11), 1413–1420 (2002) 25. Hsiao, I.T., Rangarajan, A., Gindi, G.R.: A new convex edgepreserving median prior with applications to tomography. IEEE Trans. Med. Imaging 22(5), 580–585 (2003)

146


26. Fessler, J.A.: Aspire 3.0 user’s guide: a sparse reconstruction library. Communication & Signal Processing Laboratory Technical Report No. 293, Department of Electrical and Computer Engineering, University of Michigan, Ann Arbor (1998) 27. Chang, J.-H., Anderson, J.M.M.: Regularized image reconstruction algorithms for positron emission tomography. IEEE Trans. Med. Imaging 23(9), 1165–1175 (2004) 28. Houdson, H.M., Larkin, R.S.: Accelerated image reconstruction using ordered sunsets of projection data. IEEE Trans. Med. Imaging 13(4), 601–609 (1994) 29. Fessler, J.A., Hero, A.O.: Penalized maximum-likelihood image reconstruction using space-alternating generalized EM algorithms. IEEE Trans. Image Process. 4, 1417–1429 (1995)

Yang Chen received his B.S., M.S. and Ph.D. degrees in biomedical engineering from the Department of Biomedical Engineering, First Military Medical University, china. He is currently a Lecturer at Laboratory of Image Science and Technology, Southeast University, Nanjing, China. His research interests include medical imaging reconstruction (especially on PET/CT), and medical image analysis.

Jianhua Ma received the B.S. degree in applied mathematics from the Qu-Fu Normal University, China, in 2002, the M.S. degree in computational mathematics from Sun Yat-Sen University, China, in 2005. He is currently pursuing the Ph.D. degree at the School of Biomedical Engineering, Southern Medical University. His research interests are in computed tomography imaging and inverse theory.

Qianjin Feng received his B.S., M.S. and Ph.D. degrees in biomedical engineering from the Department of Biomedical Engineering, First Military Medical University, china. He is currently an Associate Professor at Southern Medical University. His research interests include medical imaging instrumentation, image reconstruction, and medical imaging analysis.

Limin Luo is currently an professor at Laboratory of Image Science and Technology, Southeast University, Nanjing, China. He is also a senior member of IEEE. His research interests include medical imaging instrumentation and biomedical engineering.

Wufan Chen received the B.S. and M.S. degrees in applied mathematics, computational fluid dynamics from Peking University of Aeronautics and Astronautics (BUAA), China, in 1975 and 1981, respectively. From May 1981 to March 1987, he was with the School of Aerospace, National University of Defense Technology, China. From May 1987 to August 2004, he was with the Department of Training, First Military Medical University, China. Since September 2004, he has been with Southern Medical University, China, where he holds the rank of Professor in the School of Biomedical Engineering and the director of the Key Lab for Medical Image Processing of Guangdong province. His research focuses on the medical imaging and medical image analysis.