Improving image quality of x-ray in-line phase contrast imaging using an image restoration method

Xue-jun Guo; Xiao-lin Liu; Chen Ni; Bo Liu; Shi-ming Huang; Mu Gu

doi:10.1364/OE.19.023460

1. Introduction

For its great potential of enhancing x-ray imaging sensitivity and tremendously reducing radiation dose to the sample, phase contrast imaging is a powerful nondestructive detection technology in medicine, biology, and materials science [1–4]. Among all phase contrast imaging techniques, x-ray in-line phase contrast imaging [5] is the simplest, since no additional complex optical components are needed. In this technique, the phase-shifted x-ray diffracts and generates dark-bright interference fringes at boundaries and interfaces of the sample. These fringes carry phase information, and it is essential for phase retrieval and quantitative imaging [6]. However, fringe images are degraded by a practical imaging system, which impedes precise reconstruction and detail identification.

The degradation arises from partial coherence of the source, finite resolution of the detector, and unavoidable noise of the system. To obtain a higher quality image, great efforts have been made to improve source size and the imaging detector’s resolution and sensitivity [7–10]. Nevertheless, it is limited by high cost and physical limitations of the acquisition hardware. An alternative method is image restoration. Suffering from oversmoothing of edges [11–13], ringing artifacts [14–16], or being computationally intensive and complex [17–20], the classical algorithms (including frequency and spatial domain algorithms) are unable to be directly applied in phase-contrast image restoration. It is reported that the restoration algorithms extended from spatial domain to wavelet domain can deal with the oversmoothing of edge and ringing artifacts, since wavelet analysis has the properties of multiresolution and excellent time-frequency localization [21–24]. However, wavelet algorithms are still complex or slow in computation. Therefore, a special image restoration method for in-line phase contrast imaging is needed. We noted that the MTF (modulation transform function) compensation algorithm is a fast and simple tool for image restoration, although it will lead to ringing artifacts and noise amplification. If we extend the MTF compensation algorithm from Fourier domain to wavelet domain, it may provide a unique solution to phase-contrast image restoration.

In this paper, a modified algorithm combining MTF compensation with wavelet transform is proposed. Its procedures are summarized in the following. First, a series of denoising processing that preserve edges is implemented to obtain the image with high signal-to-noise ratio (SNR). Then the image is deblurred through MTF compensation in wavelet domain. Finally, the image denoising with wavelet thresholding is applied again to suppress the noise that is amplified by deblurring processing. The results indicate that the lost detail and edge information are restored effectively and artifacts are suppressed as well.

2. Degradation model of x-ray in-line phase contrast image and MTF compensation

To expound the effect of the degradation on an ideal image, the degradation model is described in a way that differs from an ideal imaging system to a practical one. The simplified schematic of x-ray in-line phase contrast imaging is shown in Fig. 1 . Consider a sample of light element material irradiated by a monochromatic x-ray of wavelength $λ$ . In the ideal imaging system, the incident wave emitted from a point source interacts with an object forming a particular distribution of wave amplitude just after the object plane. Then the wave propagates through free space and is registered in a detector. The propagation process is accompanied by the interference of scattered and transmitted waves causing the fringe along the edge of fine structure in the sample. Here we assume that the detector is an ideal detector with infinite resolution. The Fourier transform of intensity distribution in the detector plane $F (u, v)$ is calculated as follows [3]:

F (u, v) = \frac{I_{0}}{R^{2}_{1}} {1 - \sin [π λ \frac{R_{1} R_{2}}{R_{1} + R_{2}} (u^{2} + v^{2})] Φ (\frac{u}{M}, \frac{v}{M})},

where R₁ and R₂ are the distance between the source and sample, and the sample-to-detector distance, I₀; M = (R₁ + R₂)/R₁ and

Φ (\frac{u}{M}, \frac{v}{M})

denote intensity distribution of the x-ray in the source plane, magnification, and the distribution of phase change in the object plane, respectively; and u and v are spatial frequencies.

Fig. 1 Schematic of the in-line phase contrast imaging system includes an x-ray source, sample, and a detector. The profile of intensity distribution is also shown.

Download Full Size | PDF

In practice, an ideal interference pattern is always deteriorated owing to the corruption of blurring and noise. The blurring effect is caused by partial coherence of the source and spatial resolution of the x-ray detector. It can be described with MTF of the system, which is generally measured by blade-edge method [25]. The noise includes additive and multiplicative noises. The additive noise is composed of dark current and thermal noises, while the multiplicative noise is arising from the nonuniform distribution of x-ray and difference in gain coefficient of pixels.

As a partial coherence imaging system, x-ray in-line phase contrast imaging is assumed to be a linear shift-invariant system for objects with weak absorption. In frequency domain, the degradation model, i.e., the relationship among observed image $G (u, v),$ original image $F (u, v),$ additive noise image $N_{a} (u, v),$ and multiplicative noise image $N_{m} (u, v),$ is expressed as [23]

G (u, v) = k \cdot F (u, v) \cdot MTF (u, v) * N_{m} (u, v) + N_{a} (u, v),

where k and the asterisk denote the normalization factor and convolution, respectively. All noises in the system are classified as either nonstochastic or stochastic noises. The image of additive nonstochastic noise is the average of the dark current images. The image of multiplicative nonstochastic noise is obtained by

〈 {\bar{I}}_{b k} - {\bar{I}}_{d a r k} 〉

, where

G (u, v) - N (u, v)

denotes the normalization operation and

{\bar{I}}_{b k}

and

{\bar{I}}_{d a r k}

are the averages of background and dark current images, respectively.

On the basis of the degradation model, the lost information is recovered via denoising and deblurring. The additive and multiplicative nonstochastic noises are canceled respectively by subtraction and division in spatial domain. Using a logarithm, the stochastic multiplicative noise is converted to stochastic additive noise, which is erased by utilizing wavelet threshold denoise [26]. These denoising procedures effectively preserve the edge information and tremendously improve the SNR. Then a fast deblurring algorithm, which is valid only for an image with a high SNR is employed. The basic principle of the algorithm is to measure the descend degree of the MTF of the system, then calculate the ascend degree of high frequency part of the true image. To avoid noise amplification, wavelet denoise is employed again for the deblurred image. In conclusion, the restored image using MTF compensation in Fourier domain is given by [27]

I^{'} (u, v) = \frac{G (u, v) - N (u, v)}{MTF (u, v)},

where

G (u, v) - N (u, v)

denotes denoising processing. As a frequency algorithm, it is fast and precise for restoring the lost information from degenerated image, but it cannot be directly applied in x-ray in-line phase contrast imaging because of ringing artifacts and noise amplification. Thus, we extend this method from Fourier domain to wavelet domain to develop a modified algorithm for image restoration.

3. Wavelet transform theory and the modified algorithm

3.1 Wavelet transform based on FFT

Wavelet transform represents an image as a sum of wavelet functions (wavelet bases) with different locations and scales. It is implemented by a pair of analysis filters (low pass filter L and high pass filter H). Figure 2(a) shows the wavelet transform of cameraman image (256 × 256) up to the second scale, and Fig. 2(b) gives the configuration of subbands. Using the synthesis filters, the original image is exactly reconstructed from the coefficients of the bases.

Fig. 2 (a) Cameraman decomposed in seven subbands and (b) the configuration of subbands.

Download Full Size | PDF

The fast wavelet transform is conducted with the Mallat algorithm [28], which makes use of the coefficients to decompose or reconstruct the image without knowing the concrete form of scaling and wavelet functions. Although halving the computation burden, the algorithm based on the convolution or matrix vector is still computationally intensive for dealing with large size images [22]. To reduce calculation load, we can employ fast Fourier transform (FFT) to implement wavelet transform.

Let us assume that the Jth level image has a size of $\vec{i}$ , where J = 1, 2, 3... is the level of the decomposition and $\vec{j}$ is the size of origin image. The procedures of image wavelet decomposition based on FFT include the following:

(1) Consider a vector sequence $x (n) = h (n) \vec{i} + g (n) \vec{j},$

where h(n), g(n) are real coefficients with same length, corresponding to L and H filters for n = 1, 2, …

K / 2^{J}

. In Eq. (4),

\vec{i}

and

\vec{j}

denote the horizontal and vertical unit vectors, respectively. Note that L and H are orthogonal mirror filters and n × 1 matrix. The Fourier transform (FT) of Eq. (4) is expressed as

x (N) = h (N) \vec{i} + g (N) \vec{j},

where h(N) and g(N) are corresponding to h(n) and g(n) for N = 1, 2, …

K / 2^{J}

. The transposed x(N) is given by

x^{T} (N) = h^{T} (N) \vec{i} + g^{T} (N) \vec{j} .

Let

C_{0} (u, v) = F (u, v),

where C₀(u, v) denotes the FT of the finest scaling coefficients for u, v = 1, 2,…

Y_{J} (u, v) = C_{J - 1} (u, v) \cdot W (u, v)

.

(2) Two-dimensional filter coefficients W(u, v) are denoted as $W (u, v) = h (N) \otimes h^{T} (N) \vec{i} \vec{i} + h (N) \otimes g^{T} (N) \vec{i} \vec{j} + g (N) \otimes h^{T} (N) \vec{j} \vec{i} + g (N) \otimes g^{T} (N) \vec{j} \vec{j},$

where

\otimes

represents matrix multiplication.

(3) The Jth level decomposition for the image is represented as $Y_{J} (u, v) = C_{J - 1} (u, v) \cdot W (u, v) .$

Taking inverse FT of Eq. (9), we obtain

Y_{J} (m, n) = f_{L L} (m, n) \vec{i} \vec{i} + f_{L H} (m, n) \vec{i} \vec{j} + f_{H L} (m, n) \vec{j} \vec{i} + f_{H H} (m, n) \vec{j} \vec{j},

where

f_{L H} (m, n),

f_{H L} (m, n),

f_{H H} (m, n),

and

f_{L L} (m, n)

are the coefficients of horizontal, vertical, diagonal details, and approximation, respectively.

(4) Downsampling from $y_{L H} (m, n),$ $y_{H L} (m, n),$ and $y_{H H} (m, n),$ we get the horizontal, vertical, and diagonal detail subimages, respectively. The approximation subimage is obtained via downsampling from $y_{L L} (m, n)$ . Then let $C_{J} (m, n) = f_{L L} (m, n),$

Steps (1)–(4) are repeated until the required level decomposition can be obtained.

3.2 The modified algorithm for image restoration

The decomposition mentioned above can transform any single-channel linear space-invariant filtering problem into a multichannel one. Transforming Eq. (3) into wavelet yields

{\tilde{I}}_{a b}^{J} (u, v) = \frac{{\tilde{G}}_{a b}^{J} (u, v) - \tilde{N} (u, v)}{{MTF}_{L L}^{J} (u, v)},

where

G (u, v)

and

{\tilde{G}}_{a b}^{J} (u, v)

(

a, b \in L, H

) are subimages of

I (u, v)

and

G (u, v)

in wavelet domain and

{MTF}_{L L}^{J} (u, v)

is obtained from the wavelet transform of edge-blade image [27]. We get the deblurred image by inverse wavelet transform for these subimages. In Eq. (12), the assumed stationary length in this image model is reduced from

K

to

K / J

. It means that the assumption of global space-invariant is relaxed [22], leading to some “traditional” artifacts in stationary image restoration, such as ringing artifacts and oversmoothing of edges, are reduced. Certainly, too much decomposition will make the computation become complex. Thus, it is necessary to choose a proper decomposition level. The procedures of the approach are given in Fig. 3 .

Fig. 3 Flowchart of the modified MTF compensation algorithm. FWT denotes the fast wavelet transform based on FT, MTFC represents the MTF compensation algorithm, and IFWT denotes the inverse transform of FWT.

Download Full Size | PDF

4. Experiment and image quality assessment

Our modified algorithm was tested by the experimental image of a cockroach’s leg. The image acquisition system includes a micro-focus x-ray source L7901-01 of Hamamatsu photonics and an x-ray FDI cooled scientific grade detector with pixel size of 6.5 μm × 6.5 μm. The sample was located at a distance of 200 mm from the detector, the imaging magnification M was 2, and the exposure time was 30 s. The experiment was conducted with tube currents of 50 μA and tube voltage of 80 kV. Figure 4(a) illustrates the x-ray image of an edge device (1384 × 1032) made of lead, and Fig. 4(b) shows the wavelet transform of the image whose approximation subimage is adopted to compute the MTF curve. The results of the edge spread function and MTF are given in Fig. 5 .

Fig. 4 (a) X-ray image of an edge device and (b) its wavelet transform image.

Download Full Size | PDF

Fig. 5 Plots of the (a) edge spread function and (b) the MTF of the system.

Download Full Size | PDF

Figure 6(a) shows the x-ray in-line phase contrast image of the sample (1032 × 1032) acquired under the same experimental parameters. Obviously, the details and fringes nearly cannot be distinguished owing to the corruption of blur and noise. The denoising image is given in Fig. 6(b), where the details become cognizable. Utilizing the conventional MTF compensation algorithm, we obtain a deblurred image in which the lost fringe and detail information are restored. However, the deblurred image still suffers from ringing artifacts (Fig. 6(c)). From Fig. 6(d), it can be observed clearly that the lost edge and detail information are recovered without the oscillation along the edges, which demonstrates the validity of the modified MTF compensation algorithm.

Fig. 6 The results of a cockroach’s leg image: (a) original, (b) denoising, (c) conventional MTF compensation algorithm, (d) modified MTF compensation algorithm.

Download Full Size | PDF

The deblurring performance of the modified algorithm was also assessed with the parameters of image quality evaluation [29,30] by comparison between denoised and restored images, and the results are shown in Table 1 . The contrast of the image is the difference in visual properties that makes an object distinguishable from the background, and it is determined by $\sum_{i} \sum_{j} {(i - j)}^{2} \hat{p} (i, j)$ , where $\hat{p} (i, j)$ is the gray level co-occurrence matrix and i and j are the values of gray level. The contrast of restored image increases 2.17 times (from 0.06 to 0.19) compared with that of denoised image, which implies that the visibility of the image is enhanced remarkably. The energy describes the uniformity of gray level distribution within the image and is expressed as $\sum_{i} \sum_{j} {\hat{p}}^{2} (i, j)$ . Usually the image with high energy has a coarse and blurred texture. The decrease of the energy of the restored image (from 0.37 to 0.33) indicates that the texture becomes sharp. The entropy of image is a measure of information content and is defined as $\sum_{i} \sum_{j} \hat{p} (i, j) In (\hat{p} (i, j))$ . The increase of the entropy of the restored image (from 1.48 to 1.87) results in an image with richer information contents. In other words, the modified algorithm can restore the lost information effectively. The rich degree of detail signal of an image is characterized by the local energy, which is expressed as [30]

E_{l o c a l} = (\sum_{c = - 1}^{1} \sum_{d = - 1}^{1} {[f (x + c, y + d) - (\sum_{c = - 1}^{1} \sum_{d = - 1}^{1} f (x + c, y + d)) / 9]}^{2}) / {(2 K + 1)}^{2},

where

f (x, y)

denotes the gray level of an image with a size of

K^{2}

and x and y are the spatial coordinates. The variation of local energy from 0.12 to 1.11 suggests that the detail information has been recovered with our modified algorithm. The edge energy of an image represents the rich degree and clarity of edge. It is defined as

(\sum_{x} \sum_{y} e^{2} (x, y)) / K^{2}

, where

e (x, y)

is the edge image obtained by edge operator. The increase of edge energy from 0.26 to 1.97 also reveals that the lost edge information has been recovered. The variance reflects the gray level information of an image. The increase of the restored image (from 17.25 to 19.97) demonstrates that its gray level information becomes richer.

Table 1. Comparison of the Quality between Denoised and Restored Images

View Table

5. Conclusions

In conclusion, we presented a useful image restoration method for an x-ray in-line phase contrast imaging system. It extends the MTF compensation algorithm from Fourier domain to wavelet domain and is tested by an experimental image. The deblurring performance of the modified algorithm is characterized in terms of parameters of image quality evaluation. The results indicate that the contrast, edge energy, and local energy increase remarkably, which demonstrates the effectiveness of our modified algorithm. The modified algorithm will be an effective preprocessing approach for quantitative phase-contrast imaging and tomography owing to its advantages of suppressing artifacts, preserving edge information, and having low computation complexity. It can also be applied in x-ray absorption contrast imaging.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant Nos. 10875085, 91022002, and 10974143) and the Innovation Program of Shanghai Municipal Education Commission (Grant No. 11ZZ29).

References and links

1. B. Zoofan, J. Y. Kim, S. I. Rokhlin, and G. S. Frankel, “Phase-contrast x-ray imaging for nondestructive evaluation of materials,” J. Appl. Phys. 100(1), 014502 (2006). [CrossRef]

2. R. A. Lewis, “Medical phase contrast x-ray imaging: current status and future prospects,” Phys. Med. Biol. 49(16), 3573–3583 (2004). [CrossRef] [PubMed]

3. X. Wu and H. Liu, “Clinical implementation of x-ray phase-contrast imaging: theoretical foundations and design considerations,” Med. Phys. 30(8), 2169–2179 (2003). [CrossRef] [PubMed]

4. Y. S. Kashyap, P. S. Yadav, T. Roy, P. S. Sarkar, M. Shukla, and A. Sinha, “Laboratory-based X-ray phase-contrast imaging technique for material and medical science applications,” Appl. Radiat. Isot. 66(8), 1083–1090 (2008). [CrossRef] [PubMed]

5. S. Wilkins, T. Gureyev, D. Gao, A. Pogany, and A. Stevenson, “Phase contrast imaging using polychromatic hard x-ray,” Nature 384(6607), 335–338 (1996). [CrossRef]

6. X. Wu and H. Liu, “X-Ray cone-beam phase tomography formulas based on phase-attenuation duality,” Opt. Express 13(16), 6000–6014 (2005). [CrossRef] [PubMed]

7. R. Toth, J. C. Kieffer, S. Fourmaux, T. Ozaki, and A. Krol, “In-line phase-contrast imaging with a laser-based hard x-ray source,” Rev. Sci. Instrum. 76(8), 083701 (2005). [CrossRef]

8. L. Chen, L. Zheng, Y. Ai-Min, and L. Cheng-Quan, “Influence of tube voltage and current on in-line phase contrast imaging using a microfocus x-ray source,” Chin. Phys. 16(8), 2319–2324 (2007). [CrossRef]

9. Y. I. Nesterets, S. W. Wilkins, T. E. Gureyev, A. Pogany, and A. W. Stevenson, “On the optimization of experimental parameters for x-ray in-line phase-contrast imaging,” Rev. Sci. Instrum. 76(9), 093706 (2005). [CrossRef]

10. X. Wu, H. Liu, and A. Yan, “Optimization of X-ray phase-contrast imaging based on in-line holography,” Nucl. Instrum. Meth. B 234(4), 563–572 (2005). [CrossRef]

11. W. Clem Karl, “Regularization in image restoration and reconstruction,” in Handbook of Image and Video Processing, Second Edition, (Elsevier, 2005).

12. Y. T. Zhou, R. Chellappa, A. Vaid, and B. K. Jenkins, “Image restoration using a neural network,” IEEE Trans. Acoust. Speech Signal Process. 36(7), 1141–1151 (1988). [CrossRef]

13. J. K. Paik and A. K. Katsaggelos, “Image restoration using a modified Hopfield network,” IEEE Trans. Image Process. 1(1), 49–63 (1992). [CrossRef] [PubMed]

14. L. Rudin, S. Osher, and E. Fatemi, “Nonlinear total variation based noise removal algorithm,” Physica D 60(1-4), 259–268 (1992). [CrossRef]

15. T. Chan, S. Esedoglu, F. Park, and A. Yip, “Recent development of total variation image restoration,” in Handbook of Mathematical Models of Computer Vision (Springer Verlag, 2005), pp. 17–32.

16. A. Chambolle and P. Lions, “Image recovery via total variation minimization and related problems,” Numer. Math. 76(2), 167–188 (1997). [CrossRef]

17. S. Kim, “PDE-based image restoration: a hybrid model and color image denoising,” IEEE Trans. Image Process. 15(5), 1163–1170 (2006). [CrossRef] [PubMed]

18. W. Wu and A. Kundu, “Image estimation using fast modified reduced update Kalman filter,” IEEE Trans. Signal Process. 40(4), 915–926 (1992). [CrossRef]

19. V. Caselles, J. M. Morel, A. Tannenbaum, and G. Sapiro, “Introduction to the special issue on partial differential equations and geometry-driven diffusion in image processing and analysis,” IEEE Trans. Image Process. 7(3), 269–273 (1998). [CrossRef] [PubMed]

20. W. Chen, M. Chen, and J. Zhou, “Adaptively regularized constrained total least-square image restoration,” IEEE Trans. Image Process. 9(4), 589–596 (2000).

21. H. Lee and J. Paik, “Space-frequency adaptive image restoration based on wavelet decomposition,” in Proceedings of IEEE Asia Pacific Conference on Circuits and Systems (Institute of Electrical and Electronics Engineers, Seoul, 1996).

22. M. R. Banham, N. P. Galatsanos, H. L. Gonzalez, and A. K. Katsaggelos, “Multichannel restoration of single channel images using a wavelet-based subband decomposition,” IEEE Trans. Image Process. 3(6), 821–833 (1994). [CrossRef] [PubMed]

23. R. Neelamani, H. Choi, and R. Baraniuk, “ForWaRD: Fourier-wavelet regularized deconvolution for ill-conditioned systems,” IEEE Trans. Image Process. 52(2), 418–433 (2004).

24. L. Sendur and L. Selesnick, “Bivariate shrinkage functions for wavelet-based denoising exploiting interscale dependency,” IEEE Trans. Signal Process. 50(11), 2744–2756 (2002). [CrossRef]

25. E. Samei, M. J. Flynn, and D. A. Reimann, “A method for measuring the presampled MTF of digital radiographic systems using an edge test device,” Med. Phys. 25(1), 102–113 (1998). [CrossRef] [PubMed]

26. D. Donoho, “De-noising by soft thresholding,” IEEE Trans. Inf. Theory 41(3), 613–627 (1995). [CrossRef]

27. Z. Wang, Z. Geng, Y. Zhang, and X. Sui, “The MTF measurement of remote sensors and image restoration based on wavelet transform,” in Proceedings of International Conference on Wavelet Analysis and Pattern Recognition (Institute of Electrical and Electronics Engineers, Beijing, 2007), pp. 1921 −1924.

28. S. G. Mallat, “A theory for multiresolution signal decomposition: the wavelet representation,” IEEE Trans. Pattern Anal. Mach. Intell. 11(7), 674–693 (1989). [CrossRef]

29. C. H. Chen, L. F. Pau, and P. S. P. Wang, Handbook of Pattern Recognition and Computer Vision (World Science Publishing Co. Pte. Ltd, 1993), Chap 2.

30. J. Wang, X. Qi, and X. Li, “Edge and local energy NSCT based remote sensing image fusion,” J. Graduate School Chin. Acad. Sci. 26(5), 657–662 (2009) (in Chinese).

Improving image quality of x-ray in-line phase contrast imaging using an image restoration method

Abstract

1. Introduction

2. Degradation model of x-ray in-line phase contrast image and MTF compensation

3. Wavelet transform theory and the modified algorithm

3.1 Wavelet transform based on FFT

3.2 The modified algorithm for image restoration

4. Experiment and image quality assessment

5. Conclusions

Acknowledgments

References and links

Cited By

Figures (6)

Tables (1)

Equations (13)

Optics Express

	Contrast	Edge energy	Energy	Entropy	Local energy	Variance
Denoised image	0.06	0.26	0.37	1.48	0.12	17.25
Restored image	0.19	1.97	0.33	1.87	1.11	19.97