Iterative image reconstruction using non-local means with total variation from insufficient projection data

Abstract

In this work, algebraic reconstruction technique (ART) is extended by using non-local means (NLM) and total variation (TV) for reduction of artifacts that are due to insufficient projection data. TV and NLM algorithms use different image models and their application in tandem becomes a powerful denoising method that reduces erroneous variations in the image while preserving edges and details. Simulations were performed on a widely used 2D Shepp-Logan phantom to demonstrate performance of the introduced method (ART+TV) NLM and compare it to TV based ART (ART+TV) and ART. The results indicate that (ART+TV) NLM achieves better reconstructions compared to (ART+TV) and ART.

Keywords

Limited-angle tomography total variation non-local means insufficient projection

1 Introduction

X-ray dose delivered to the patient has been a big concern in computed tomography (CT) imaging. Realistically, a full abdominal CT scan gives a typical effective dose of 10.0 mSv which is 3.3 years of equivalent effective dose exposed from natural background radiation while a chest X-ray is only 2.4 days [1]. As a result of this important issue, reducing the amount of exposed radiation while keeping acceptable image resolution has become an active research area. Lowering the dose itself not only reduces the amount of patient dose but also reduces image resolution which is a crucial parameter in CT imaging. CT imaging based on limited-view angle [2 –4] and fewer projections [4, 5] have been used to reduce X-ray dose. However, as the number of projections decreases, data may become insufficient to fully recover the exact image. Analytical methods such as Fourier transform (FT) or filtered back projection (FBP) create streaking artifacts and distortions appearing in the reconstructed image.

In order to overcome this problem, iterative reconstruction techniques and image models have been utilized. Algebraic reconstruction technique (ART) has been one of the commonly used iterative method to estimate the image from measured projection data [4]. A simultaneous form of ART (SART) was also proven to have satisfactory results in tomosynthesis [5] and limited view angle imaging [6].

Besides the image reconstruction method, a prior image model is typically used to improve the inconsistency, error, and insufficient information in CT imaging. After its introduction by Rudin, Osher and Fatemi (ROF), total variation (TV) has been used commonly as an image model [7]. Due to this interest, many variations of TV such as “anisotropic” [8], “nonlocal” [9], “nonlinear” [10] TV were introduced.

In parallel to develop new image models, effective solutions for TV models were introduced [11, 12]. Following these developments TV model was used in many image processing application such as denoising [13 –15], inpainting [16 –20], deconvolution [21 –23], and image restoration [24 –29]. Following these advances, variational techniques have been adapted to medical imaging field to improve existing tomographic image reconstruction methods. TV based algorithms were especially useful for limited-angle tomography [30], dose reduction [31], and digital breast tomosynthesis [32]. Similarly it has been used for other medical imaging modalities such as magnetic resonance imaging (MRI) [33], positron emission tomography (PET) [34] etc.

Besides TV, non-local means (NLM) filtering algorithm has commonly been used in image denoising for preserving fine details [35] and iterative sparse-angular CT reconstruction [36]. The underlying principle of NLM is the similarity of nonlocal patches in the images. Hence TV and NLM use different prior information about images: Sparsity of gradients (TV), similarity of local patches (NLM). Hence, it is intuitive to design better algorithms that integrate NLM and TV into iterative reconstruction algorithm.

The contribution of this work is the usage of both TV and NLM as image prior in order to improve tomographic reconstruction. A similar approach was used for image enhancement [37]. However to the best of our knowledge, it has not been used in tomographic reconstruction yet.

In this study, a widely used iterative reconstruction algorithm, ART, is combined with TV and NLM which will be named as (ART+TV)NLM in the paper from now on. In the proposed algorithm, the image reconstructed by the widely used ART+TV algorithm is filtered by NLM sequentially. The Shepp-Logan phantom was used in the simulations to compare the performance of the methods.

The rest of the paper is organized as follows: We briefly described the conventional reconstruction methods and denoising tools in the following Section. Then experimental setup was detailed and comparisons of (ART+TV)NLM method with the other methods were shown on Shepp-Logan both qualitatively and quantitatively in Section III. Finally, the paper was concluded in Section IV.

2 Methods

Image reconstruction problem studied in this work can be modeled by using the following formulation: $Y_{i} = \sum_{j = 1}^{N} A_{ij} X_{j} i = 1, 2, \dots, M,$ (1) where, Y _i is the ray-sum computed for the ith ray while X _j represents the intensity value of the jth pixel. A _ij is the contribution factor of the pixel X _j with respect to the ray-sum Y _i, M and N refer to the total number of ray-sums and pixels respectively. Iterative methods adopt a numerical approximation approach where the measured (detector images) and calculated projection are compared and the difference is then backward projected over the related rays to update the current image until a specific convergence is satisfied. The following formulation shows how an image is updated by using the ARTalgorithm: $X_{j}^{(k + 1)} = X_{j}^{(k)} + \frac{Y_{i} - \sum_{k = 1}^{N} A_{ik} X_{k}^{(k)}}{\sum_{k = 1}^{N} A_{ik}^{2}} A_{ij}, \begin{matrix} i = 1, 2, \dots, M \\ j = 1, 2, \dots, N \end{matrix}$ (2) where, $X_{j}^{(k)}$ and $X_{j}^{(k + 1)}$ show the value of pixel intensities in the current and next iterations respectively. The algorithm Equation (2) consists of three major steps: The forward projection of the estimated image, calculation of the difference between the measured and calculated projection data, and the backward projection of the difference to the image according to the pixel contribution with respect to the related rays. This process is repeated until all projections are considered and the estimated image converges to a solution.

However, reconstructing an exact image from a highly missing projection data set by using iterative reconstruction shows distortions and streaking artifacts; such as, salt and pepper noise. TV minimization has been a solution to acquire better image quality while preserving edges and reducing the background noise. TV shows how smoothly the intensity of an image changes. In most tomographic images, the intensities are relatively constant within organs or large tissues. Hence minimizing the TV of an image reduces the inhomogeneity within organs and tissues while preserving edges. In order to adapt this into an iterative reconstruction, ℓ₁ norm of the gradient image is minimized. This method is also known as the total variation (TV) method and formulated as follows for a 2D image: $T V (X) = \sum_{i}^{K} \sum_{j}^{L} \sqrt{{(X_{i, j} - X_{i - 1, j})}^{2} + {(X_{i, j} - X_{i, j - 1})}^{2}}$ (3)

NLM has also been utilized in image processing for denoising. The denoising process is repeated pixel by pixel for the entire image and formulated as: $N L M (X_{i}) = \sum_{X_{j \in S W}} W (X_{i,} X_{j}) X_{j}$ (4) where X _i and X _j are the intensity values of image X at pixels i and j. X _j is limited to a search window (SW) which bounds the neighboring pixels remoteness. μ denotes the intensity value of specific pixel at image X. W (X _i, X _j) represents the weighting function between pixels i and j. The weighting function shows how much the difference between pixels i and j is penalized and it is calculated by using the following formulation: $w (X_{i}, X_{j}) = \sum_{δ \in P} e^{- \frac{G_{a} * {| (x_{i} + δ) - (x_{j} + δ) |}^{2}}{h^{2}}}$ (5) where, G _a is Gaussian kernel, h is the filtering parameter which controls the power of the filter, and ^* denotes convolution. Filtering parameter h is directly related to the level of noise in the image. Symbol δ represents the patch (P) centered at pixels i and j. The weights are normalized as follows: $C (X) = \sum_{X_{J} \in SW} w (X_{i}, X_{j})$ (6)

The main focus of this study is to combine two powerful noise reduction steps to achieve a higher quality in the reconstructed image. TV assumes the sparsity of the gradients whereas NLM assumes the repetition of patches within images. As both assumptions are correct and separate from each other, their combination within an algorithm will benefit from both image models. Hence, a new approach that we name (ART+TV)NLM is proposed. ART+TV has been widely accepted as a powerful iterative image reconstruction algorithm with a feature of preserving edges while reducing the background noise. However small details can be distorted by the power of TV, thus image reconstructed by ART+TV was then filtered by NLM to make smaller details more visible.

3 Simulation setup and numerical results

3.1 Simulation setup

Simulation of (ART+TV)NLM was performed on a widely used Shepp-Logan phantom with a resolution of 128×128 pixels. The detector with a size of 256 was located 164 pixels away from the center of the object. The X-ray point was located 512 pixels above the center of the object. The motion of X-ray is in step and shoot manner while detector rotates opposite to the rotation direction of X-ray to acquire an arc-like position in the simulations. The most important part of all limited view angle and fewer projection problems is the scan angle and the number of projection as they have a direct impact on the image quality. In this study the scan angle was chosen to be 180° covering the half of the full rotation with an increment of 9°. Thus the projection number is limited to 21.

The choice of the parameters which control the power of TV and NLM during the reconstruction process is critical as it affects reconstruction quality. The regularization parameter of TV which controls the power of TV inclusion was set constant to a manually adjusted value of 0.8 for all simulations. Search window and patch size for NLM filtering were selected as 15, 11 respectively as these values were optimized for 128×128 Shepp-Logan phantom in our previous work [38]. The most important parameter in NLM filtering is the filtering parameter, h, which is directly related to the noise level of the reconstructed image. Thus, appropriate selection of the filtering parameter affects the performance of the filter. We know from our previous studies that in limited-view angle iterative reconstructions, the root mean square error (RMSE) decreases exponentially with the number of iterations [39, 40]. The decrease in RMSE is sharper in the first 4-5th iterations while the slope of the curve slightly decreases to a stable form until the 15th iteration with almost no further important change occurs after the 15th iteration. Considering these details, the filtering parameter h was adaptively changed according to iteration number in an exponential decay form rather than using a stable h value. The effect of choosing an exponential h instead of constant h is discussed in the following part.

3.2 Visual and numerical assessments

Performance analysis of (ART+TV)NLM method was assessed both qualitatively by visual assessment and quantitatively by using RMSE, SNR and SSIM parameters. The reconstruction results on Shepp-Logan phantom are presented in Fig. 1. The figure shows the original phantom and images reconstructed with FBP, ART, ART+TV and (ART+TV)NLM respectively. It can be clearly seen that inclusion of NLM increased the image quality and the proposed method has better results than ART+TV. (ART+TV)NLM algorithm performs a very close result with the original phantom with almost negligible artifacts in the image. To further visualize the difference between reconstruction algorithms, two different region of interests (ROIs) were selected and zoomed in to clearly show the superior performance of the proposed method over the other methods. In Fig. 2, FBP and ART cannot yield satisfactory results as the artifact level makes small details indistinguishable. However the inclusion of TV makes small details more visible while still having undesired artifacts at the boundaries. Combining NLM with TV gives the closest match with the original image by removing the existent background noise and undesired artifacts in ART+TV.

The appropriate selection of the filtering parameter should be taken into account as it directly affects the performance of the reconstruction. In order to show this impact, the effect of the filtering parameter h on the reconstructed image was evaluated by selecting the parameter constant and exponential decay as described above and their results are shown in Fig. 3. The image reconstructed with a constant h showed lots of artifacts by generating blurs and removing the fine details out whereas exponential h helped in obtaining a much better reconstructed image as shown in Fig. 3

To visualize the difference between ART+TV and (ART+TV)NLM, both vertical and horizontal profiles were drawn from the reconstructed images. The vertical profile was from 20th row to the 107th row at 64th column while the horizontal profile was extracted from 40th to 90th column at 102th row as shown in Fig. 4. The original phantom profiles were added at the same selected rows and columns as a reference for the reconstruction results. The profile comparison results with the reference profiles clearly show that the inclusion of NLM provides remarkable improvement over ART+TV especially at the boundaries where the total variation is higher. With its superior background noise removal feature, the NLM perfectly approaches to the original profiles while ART+TV results show fluctuations with sharp spikes at the uniform areas.

For a deeper performance analysis RMSE values for the reconstructions are presented in Fig. 5. (ART+TV) NLM gives slightly the lowest RMSE value among three reconstruction methods. In order to increase the consistency with RMSE results, one of the well-known quality assessment metric, Structure SIMilarity (SSIM) index was utilized. SSIM index is utilized to mimic the numerical analysis for human vision system. SSIM parameter varies from 0 to 1 as it gives the best image at value “1”. In Fig. 5, (ART+TV)NLM achieved the exact value shows the full reconstruction “1” while ART and ART+TV reached to a value of 0.774 and 0.989 at the 20th iteration.

The last parameter used for the numerical comparison is the SNR value. Figure 6 gives the SNR graph for the reconstruction algorithms. (ART+TV)NLM method achieved the highest SNR with a value of 45.81 while ART and ART+TV reached to 19.24 and 33.21 respectively at the 20th iteration.

4 Conclusion

In this paper, an iterative reconstruction algorithm based on combining TV minimization and NLM filtering has been proposed for insufficient projection data problem. The numerical results were conducted to compare the performances of ART, ART+TV, (ART+TV)NLM by using a widely used Shepp-Logan phantom. The introduced algorithm, (ART+TV)NLM, showed better results than the other two reconstruction methods in terms of RMSE, SSIM, and SNR metrics.

Footnotes

Acknowledgments

This work has been supported by TUBITAK, The Scientific and Research Council of Turkey, under the grant 111E086.

References

European Commission, Radiation Protection Report 118, “Referral guidelines for imaging”. Directorate General for the Environment of the European Commission; 2000.

Oskoui

and Stark

, A comparative study of three reconstruction methods for a limited view computer tomography problem, IEEE Trans Med Img 8(1) (1989), 43–49.

Andersen

A.H.

, Algebraic reconstruction in CT from limited views, IEEE Trans Med Img 8(1) (1989), 50–55.

Kak

A.C.

and Slaney

, Principles of computerized tomographic imaging. IEEE, New York, 1988.

Chen

G.H.

, Tang

and Leng

, Prior Image constrained compressed sensing (PICCS): A method to accurately reconstruct dynamic CT images from highly undersampled projection data sets, Med Phys 35(2) (2008), 660–663.

, Khamene

and Fluck

, High performance tomosynthesis enabled via a GPU-based iterative reconstruction framework, Proc of SPIE 7258 (2009), 72585A.

Rudin

Leonid I.

, Osher

and Fatemi

, Nonlinear total variation based noise removal algorithms, Physica D: Nonlinear Phenomena 60(1) (1992), 259–268.

Selim

Esedoglu

and Osher

S.J.

, Decomposition of images by the anisotropic Rudin-Osher-Fatemi model, Communications on Pure and Applied Mathematics 57(12) (2004), 1609–1626.

Antoni

Buades

Coll

and Morel

J.-M.

, A non-local algorithm for image denoising, Computer Vision and Pattern Recognition, 2005. CVPR 2005. IEEE Computer Society Conference on, Vol. 2. IEEE, 2005.

10.

Rudin Leonid

, Osher

and Fatemi

, Nonlinear total variation based noise removal algorithms, Physica D: Nonlinear Phenomena 60(1) (1992), 259–268.

11.

Chambolle

and Lions

P.L.

, Image recovery via total variation minimization and related problems, Numerische Mathematik 76(2) (1997), 167–188.

12.

Chambolle

, An algorithm for total variation minimization and applications, Journal of Mathematical Imaging and Vision 20(1-2) (2004), 89–97.

13.

Vogel

C.R.

and Oman

M.E.

, Iterative methods for total variation denoising, SIAM Journal on Scientific Computing 17(1) (1996), 227–238.

14.

Dobson

D.C.

and Vogel

C.R.

, Convergence of an iterative method for total variation denoising, SIAM Journal on Numerical Analysis 34(5) (1997), 1779–1791.

15.

Beck

and Teboulle

, Fast gradient-based algorithms for constrained total variation image denoising and deblurring problems, Image Processing, IEEE Transactions on 18(11) (2009), 2419–2434.

16.

Chan

T.F.

, Shen

and Zhou

H.M.

, Total variation wavelet inpainting, Journal of Mathematical Imaging and Vision 25(1) (2006), 107–125.

17.

Zhang

and Chan

T.F.

, Wavelet inpainting by nonlocal total variation, Inverse Problems and Imaging 4(1) (2010), 191–210.

18.

Wen

Y.W.

, Chan

R.H.

and Yip

A.M.

, A primal–dual method for total-variation-based wavelet domain inpainting, Image Processing, IEEE Transactions on 21(1) (2012), 106–114.

19.

Bertalmio

, Vese

, Sapiro

and Osher

, Simultaneous structure and texture image inpainting, Image Processing, IEEE Transactions on 12(8) (2003), 882–889.

20.

Bertalmio

, Sapiro

, Caselles

and Ballester

, In Proceedings of the 27th Annual Conference on Computer Graphics and Interactive Techniques, ACM Press/Addison-Wesley Publishing Co., 2000, pp. 417–424.

21.

Chan

T.F.

, Yip

A.M.

and Park

F.E.

, Simultaneous total variation image inpainting and blind deconvolution, International Journal of Imaging Systems and Technology 15(1) (2005), 92–102.

22.

Bioucas-Dias

J.M.

, Figueiredo

M.A.

and Oliveira

J.P.

, Total variation-based image deconvolution: A majorization-minimization approach, In Acoustics, Speech and Signal Processing, 2006 ICASSP 2006 Proceedings 2006 IEEE International Conference on, Vol. 2, 2006, pp. II–II. IEEE.

23.

Oliveira

J.P.

, Bioucas-Dias

J.M.

and Figueiredo

M.A.

, Adaptive total variation image deblurring: A majorization–minimization approach, Signal Processing 89(9) (2009), 1683–1693.

24.

Osher

, Burger

, Goldfarb

, Xu

and Yin

, An iterative regularization method for total variation-based image restoration, Multiscale Modeling & Simulation 4(2) (2005), 460–489.

25.

Wang

, Yang

, Yin

and Zhang

, A new alternating minimization algorithm for total variation image reconstruction, SIAM Journal on Imaging Sciences 1(3) (2008), 248–272.

26.

Chan

T.F.

, Golub

G.H.

and Mulet

, A nonlinear primal-dual method for total variation-based image restoration, SIAM Journal on Scientific Computing 20(6) (1999), 1964–1977.

27.

Chan

, Esedoglu

, Park

and Yip

, Recent developments in total variation image restoration, 17, Mathematical Models of Computer Vision, 2005.

28.

Osher

, Solé

and Vese

, Image decomposition and restoration using total variation minimization and the H 1, Multiscale Modeling & Simulation 1(3) (2003), 349–370.

29.

M.K.

, Weiss

and Yuan

, Solving constrained total-variation image restoration and reconstruction problems via alternating direction methods, SIAM Journal on Scientific Computing 32(5) (2010), 2710–2736.

30.

Chen

, Jin

, Li

and Wang

, A limited-angle CT reconstruction method based on anisotropic TV minimization Phys, Med Biol 58, (2013), 2119–2141.

31.

Chen

G.H.

, Tang

and Leng

, Prior Image constrained compressed sensing (PICCS): A method to accurately reconstruct dynamic CT images from highly undersampled projection data sets, Med Phys 35(2) (2008), 660–663.

32.

Sidky

E.Y.

, Pan

, Reiser

I.S.

, Nishikawa

R.M.

, Moore

R.H.

and Kopans

D.B.

, Enhanced imaging of microcalcifications in digital breast tomosynthesis through improved image-reconstruction algorithms, Med Phys 36(11) (2009), 4920–4932.

33.

Block

K.T.

, Uecker

and Frahm

, Undersampled radial MRI with multiple coils. Iterative image reconstruction using a total variation constraint, Magnetic Resonance in Medicine 57(6) (2007), 1086–1098.

34.

Sawatzky

, et al. Accurate EM-TV algorithm in PET with low SNR, Nuclear Science Symposium Conference Record, 2008, NSS’08 IEEE IEEE, 2008.

35.

Lou

, Zhang

and Osher

, Image recovery via nonlocal operators, J Sci Compt 42(2) (2010), 185–197.

36.

Huang

, Ma

, Liu

, Zhang

, Bian

, Feng

and Chen

, Sparse angular CT reconstruction using non-local means based iterative-correction POCS, Comp Bio Med 41(4) (2011), 195–205.

37.

Likforman-Sulema

, Darbonb

and Barney

E.H.

, Smith, Enhancement of historical printed document images by combining Total Variation regularization and Non-local Means filtering, Im and Com Vis 29(5) (2011), 351–363.

38.

Ertas

, Akan

, Yildirim

and Kamasak

, Image Denoising by using Non-Local Means and Total Variation, Proc of Signal Processing and Application Conference, Trabzon, 2014.

39.

Ertas

, Yildirim

, Kamasak

and Akan

, Digital breast tomosynthesis image reconstruction using 2D and 3D total variation minimization Biomed, Eng Onl 12, (2013), 112.

40.

Sechopoulos

, A review of breast tomosynthesis. Part I. The image acquisition process, Med Phys 40(1) 2013, 014301–1:12.