Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2705164

Formats

Article sections

Authors

Related links

Opt Commun. Author manuscript; available in PMC 2010 August 15.

Published in final edited form as:

Opt Commun. 2009 August 15; 282(16): 3392–3396.

doi: 10.1016/j.optcom.2009.05.027PMCID: PMC2705164

NIHMSID: NIHMS118119

Fanbo Meng, Center for Bioengineering and School of Electrical and Computer Engineering, University of Oklahoma, Norman, OK 73019, USA;

Xizeng Wu: ude.cmbau@uwx; Hong Liu: ude.uo@uil

Previous studies have shown that iterative in-line x-ray phase retrieval algorithms may have higher precision than direct retrieval algorithms. This communication compares three iterative phase retrieval algorithms in terms of accuracy and efficiency using computer simulations. We found the Fourier transformation based algorithm (FT) is of the fastest convergence, while the Poisson-solver based algorithm (PS) has higher precision. The traditional Gerchberg-Saxton algorithm (GS) is very slow and sometimes does not converge in our tests. Then a mixed FT-PS algorithm is presented to achieve both high efficiency and high accuracy. The mixed algorithm is tested using simulated images with different noise level and experimentally obtained images of a piece of chicken breast muscle.

In the booming field of in-line phase-sensitive x-ray imaging, there are two closely related techniques: phase contrast imaging and phase retrieval [1]. The former is a imaging modality that aims to enhance the total contrast of the images using the differential phase shift induced by the object under investigation. By presenting much higher visibility of details inside soft tissues than the traditional attenuation based imaging, it is promising to be a valuable diagnostic tool for applications like mammography. The latter is basically a mathematical technique of solving the diffraction integral, which describes the formation of the phase contrast images, for the quantitative phase shift map. The phase map may not be able to give us a clearer observation of the structure of the object, but its quantitative nature can lead to essential information on the internal structures complementary to attenuation.

Researchers aimed to develop deterministic non-iterative phase retrieval algorithms by “linearize” their theoretical formulations with respect to the phase term. Recent studies [2, 3], however, showed that iterative algorithms based on a more general “nonlinear” formula can improve the precision of the resultant phase maps. In this paper, we implemented and compared three different iterative phase retrieval algorithms. First is the Fourier transformation based algorithm (abbreviated as FT algorithm hereafter), which solves for the phase term in the spatial frequency domain, second is the Poisson-solver based algorithm (PS), which solves for the phase term in the spatial domain using readily available Poisson equation solving algorithms and the last is a variant of the traditional Gerchberg-Saxton algorithm (GS) [4]. The three algorithms were compared in terms of precision and efficiency. The results of the comparison suggest us to develop a mixed algorithm with the both merits of high accuracy and computational speed. Computer simulations were utilized to investigate the robustness of the mixed algorithm against image noise. The algorithm was then tested using experimentally obtained images of a piece of chicken breast muscle.

In the physical picture, the formation of phase contrast in the in-line phase contrast imaging scheme is actually the process of x-ray wave diffraction and, therefore, can be described by wave diffraction formulations that have been intensively studied in optics. The most general formula for wave diffraction under the scalar wave approximation is the Rayleigh-Sommerfeld diffraction integral [5]

$$U(x,y;z)=\frac{1}{2\pi}\iint U({x}^{\prime},{y}^{\prime};0)\times \frac{\text{exp}(-\mathit{\text{ikr}})}{r}\frac{z}{r}\phantom{\rule{0.2em}{0ex}}\left(\mathit{\text{ik}}+\frac{1}{r}\right)\phantom{\rule{0.2em}{0ex}}d{x}^{\prime}d{y}^{\prime},$$

(1)

where *U*(*x, y; z*) is the complex valued wavefield distribution on the *z* plane, *k* is the wave number, *r* = [(*x* − *x′*)^{2} + (*y* − *y′*)^{2} + *z ^{2}*]

$$U({u}_{x},{u}_{y};z)=U({u}_{x},{u}_{y};0)\cdot \text{exp}\phantom{\rule{0.2em}{0ex}}\left[-i2\pi z\sqrt{\frac{1}{\lambda}-{u}_{x}^{2}-{u}_{y}^{2}}\right]\phantom{\rule{0.2em}{0ex}},$$

(2)

where (*u _{x}*,

When the wavefield can be approximated by paraxial waves, Equation (1) is reduced to the Fresnel-Kirchhoff integral [7]

$$\begin{array}{cc}U(x,y;z)=\frac{i\phantom{\rule{0.1em}{0ex}}\text{exp}\phantom{\rule{0.1em}{0ex}}(-\mathit{\text{ikz}})}{\lambda z}\iint \hfill & U({x}^{\prime},{y}^{\prime};0)\times \hfill \\ \hfill & \text{exp}\phantom{\rule{0.3em}{0ex}}\left[-\mathit{\text{ik}}\frac{{(x-{x}^{\prime})}^{2}+{(y-{y}^{\prime})}^{2}}{2z}\right]\phantom{\rule{0.3em}{0ex}}d{x}^{\prime}d{y}^{\prime}\hfill \end{array}$$

(3)

and Equation (2) to

$$U({u}_{x},{u}_{y};z)=U({u}_{x},{u}_{y};0)\cdot \text{exp}\phantom{\rule{0.2em}{0ex}}\left[-\mathit{\text{ikz}}-i\lambda z({u}_{x}^{2}+{u}_{y}^{2})\right].$$

(4)

If the object under investigation is placed on the plane *z* = 0 (referred to as the object plane hereafter), the wavefield on the object plane is described phenomenologically as the product of the incident wave *U _{in}*(

$$U(x,y;0)={U}_{\mathit{\text{in}}}(x,y;0)\cdot T(x,y).$$

(5)

The experimentally obtained images at a z-plane is the intensity distribution of the wavefield, i.e., the squared modulus of *U*(*x, y; z*),

$$I(x,y;z)={\left|U(x,y;z)\right|}^{2}=U(x,y;z)\cdot {U}^{\ast}(x,y;z).$$

(6)

The problem of phase retrieval is then defined in the literature as solving Equation (6) for *T*(*x*, *y*) based on two or more experimentally obtained images and a given incident wavefield *U _{in}*(

A general iterative approach for phase retrieval algorithms was first presented by Gerchberg and Saxton [4] and then further developed by other researchers [8, 9]. Now it has evolved into a large bunch of related iterative algorithms with a broad spectrum of applications in electron microscopy, wave front sensing, astronomy, crystallography and many other fields. A recent paper [10] summarized several common variant GS algorithms. Reported in the literature [11, 12, 13] were also applications of GS algorithms in the in-line x-ray phase retrieval.

The GS algorithm, utilizing the fact that wavefield propagation is reversible, projects the wavefield back and forth between the object plane and the image plane repeatedly and applies certain constraints in each iteration. In our algorithm, we take the subsidiary constraint which replace the modulus of the wavefield with the square root of the intensity measurement while leaving the phase of the wavefield unchanged. The algorithm begins with an initial guess *U*^{(0)}(*x*, *y*; 0) and then follows the iteration cycle below:

- let ${U}^{(n)}(x,y;0)\Leftarrow \sqrt{I(x,y;0)}\text{exp}\phantom{\rule{0.2em}{0ex}}\left[i\text{arg}{U}^{(n)}(x,y;0)\right];$
- let ${U}^{(n)}(x,y;{R}_{2})\Leftarrow \sqrt{I(x,y;{R}_{2})}\text{exp}\phantom{\rule{0.2em}{0ex}}\left[i\text{arg}{U}^{(n)}(x,y;{R}_{2})\right];$

Here “arg” means the argument of an complex number and superscript (*n*) and (*n*+1) means iteration count. The iteration stops when, for example, the change in *U*(*x*, *y*; 0) falls below an predefined tolerance level. From *U*(*x*, *y*; 0) and *U _{in}*(

It is worth noting that, to take into account the effect of the partial spatial coherence of the incident x-ray wavefield and the frequency response of the detector, the obtained images should be de-convoluted according to the total optical transfer function (OTF) first. Please see the next section for details about the total OTF.

A general formulation for phase retrieval was first presented by Liu and Wu [14, 15], and was then re-derived and experimentally proved in other groups' publications [16, 3]. It was shown that this formulation covers other approaches based on the transport of intensity equation (TIE) and the contrast transfer function (CTF). This formulation is based on the Fresnel-Kirchhoff approximation and was derived using the Wigner distribution [15] of the wavefield to incorporate the partial spatial coherence of *U _{in}*(

$$\begin{array}{c}\widehat{\mathcal{F}}\phantom{\rule{0.2em}{0ex}}\left[{M}^{2}I(M\overrightarrow{\rho};{R}_{2})\right]={I}_{0}\cdot \text{OTF}\phantom{\rule{0.2em}{0ex}}\left(\frac{\overrightarrow{u}}{M}\right)\phantom{\rule{0.2em}{0ex}}\{\text{cos}\phantom{\rule{0.1em}{0ex}}(\alpha {u}^{2})\widehat{\mathcal{F}}\phantom{\rule{0.2em}{0ex}}\left[{A}^{2}(\overrightarrow{\rho})\right]\hfill \\ \phantom{\rule{14em}{0ex}}+2\phantom{\rule{0.1em}{0ex}}\text{sin}\phantom{\rule{0.1em}{0ex}}\left(\alpha {u}^{2}\right)\phantom{\rule{0.1em}{0ex}}\widehat{\mathcal{F}}\phantom{\rule{0.1em}{0ex}}\left[{A}^{2}(\overrightarrow{\rho})\phantom{\rule{0.1em}{0ex}}\phi \phantom{\rule{0.1em}{0ex}}(\overrightarrow{\rho})\right]\hfill \\ \phantom{\rule{12em}{0ex}}+i\frac{\lambda {R}_{2}}{2M}sin\phantom{\rule{0.1em}{0ex}}(\alpha {u}^{2})\phantom{\rule{0.1em}{0ex}}\overrightarrow{u}\cdot \widehat{\mathcal{F}}\phantom{\rule{0.1em}{0ex}}\left[\nabla {A}^{2}(\overrightarrow{\rho})\right]\hfill \\ \phantom{\rule{12em}{0ex}}+i\frac{\lambda {R}_{2}}{M}\phantom{\rule{0.1em}{0ex}}\text{cos}\phantom{\rule{0.1em}{0ex}}(\alpha {u}^{2})\phantom{\rule{0.1em}{0ex}}\overrightarrow{u}\cdot \widehat{\mathcal{F}}\phantom{\rule{0.1em}{0ex}}\left[\phi \phantom{\rule{0.1em}{0ex}}(\overrightarrow{\rho})\nabla {A}^{2}(\overrightarrow{\rho})\right]\},\hfill \end{array}$$

(7)

where = (*x*, *y*) and *$\stackrel{\u20d7}{u}$* = (*u _{x}*,

$$\text{OTF}\left(\frac{\overrightarrow{u}}{M}\right)={\text{OTF}}_{\text{G}.\text{U}.}\phantom{\rule{0.1em}{0ex}}\left(\frac{\overrightarrow{u}}{M}\right)\cdot {\text{OTF}}_{\text{det}}\phantom{\rule{0.1em}{0ex}}\left(\frac{\overrightarrow{u}}{M}\right)$$

(8)

is the product of the OTFs of the x-ray source and the detector.

An iterative algorithm was presented based on Equation (7) [2]. This algorithm solves *ϕ*() from Equation (7) in the spatial frequency domain utilizing the fast Fourier transformation (FFT) technique. The algorithm can be described as

$$\begin{array}{c}\phantom{\rule{2em}{0ex}}\widehat{\mathcal{\text{F}}}\phantom{\rule{0.2em}{0ex}}\left[{A}^{2}\right]=\stackrel{\sim}{I}(\overrightarrow{\rho};0)/{I}_{0}\hfill \\ \widehat{\mathcal{F}}\phantom{\rule{0.2em}{0ex}}\left[{A}^{2}{\phi}^{(n+1)}\right]=\frac{{M}^{2}\stackrel{\sim}{I}(M\overrightarrow{\rho};{R}_{2})/{I}_{0}-\text{cos}(\alpha {u}^{2})\widehat{\mathcal{F}}\phantom{\rule{0.2em}{0ex}}\left[{A}^{2}\right]-T4({A}^{2},{\phi}^{(n)},{R}_{2})}{2\phantom{\rule{0.1em}{0ex}}\text{sin}\phantom{\rule{0.1em}{0ex}}(\alpha {u}^{2})},\hfill \end{array}$$

(9)

where *I˜*(*; z*) is the image de-convoluted by OTF(*$\stackrel{\u20d7}{u}$/M*), and *T*4 is the fourth term in the bracket in Equation (7).

It was shown that, when *α* 1 holds, Equation (7) can be reduced to a Poisson-type equation in the spatial domain as [14]

$${M}^{2}\stackrel{\sim}{I}(M\overrightarrow{\rho};{R}_{2})-\stackrel{\sim}{I}(\overrightarrow{\rho};0)=\frac{\lambda {R}_{2}}{2\pi M}\nabla \cdot \left[\stackrel{\sim}{I}(\overrightarrow{\rho};0)\nabla \phi \phantom{\rule{0.1em}{0ex}}(\overrightarrow{\rho})\right].$$

(10)

Equation (10) can also be obtained as an approximation of the transport of intensity equation (TIE) [17, 18] in the case of a point x-ray source.

Equation (10) can be solved using a large bunch of mature numerical algorithms generally called “Poisson solvers”, which are readily available either commercially or for free in various programming languages.

We tested the three iterative algorithms, the Fourier transformation based algorithm (FT), the Poisson-solver based algorithm (PS) and the Gerchberg-Saxton algorithm (GS), first using images produced by a computer simulation program, then images obtained in an experiment.

The three phase retrieval algorithms described above, together with a simulation algorithm that produces phase contrast images based a given *T*(*x*, *y*) distribution, were coded in c++ and Fortran programming languages and compiled into a single program. The compiling and running environment is a desktop personal computer with an Intel Pentium 4 2800 MHz CPU and 2G physical memory running the Microsoft Windows XP operating system. The FFTW library (http://www.fftw.org/) was used for fast Fourier transformation calculations. Line successive over relaxation (LSOR) [19] and alternating direction implicit (ADI) [20] methods were employed for solving the Poisson-type equation.

We cropped from a clinical mammographic image two maps (shown in Figure 1) that simulate the modulus *A*(*x*, *y*) and argument *ϕ*(*x*, *y*) of the optical transmission function of an virtual object. This virtual object resembles a piece of soft tissue typical of mammographic studies. Then the virtual object was used for computer simulations. An attenuation based image and a phase contrast image of the virtual object were produced according to Equation (2).

The attenuation and phase maps used in the computer simulations. Following the convention in x-ray radiography, in the images hereafter, black means larger value and white means smaller value.

Both maps are 500 by 500 in size. The value of map *A*(*x*, *y*) ranges from 0.802 to 0.899 with a standard deviation (SD) of 0.025; that of *ϕ*(*x*, *y*) ranges from −5.46 to 0 with a SD of 1.578.

The parameters for the simulation include: *R*_{2} = 10m, *M* = 10, photon energy= 20keV, pixel pitch of the detector is 20*μ*m. The simulated phase-contrast image *I*_{2} under ideal conditions are shown in Figure 2.

We first tested the precision and efficiency of the three phase retrieval algorithms using the simulated images. We measure the precision of the retrieval by the SD of the difference image between the retrieved and the set phase map (abbreviated as SDDI hereafter). The SDDI is plotted in Figure 3 against program running time.

We found that the FT algorithm is about two orders of magnitude faster than the PS algorithm, but the PS algorithm can result in a smaller SDDI than the FT algorithm. The GS algorithm is much slower than the FT and PS algorithms. This indicted us to try mixed algorithms by applying the FT, then the PS and finally the GS algorithm in turn in order to achieve the goal of both speed and accuracy. As can be seen from Figure 3, the PS algorithm did help improve the result of FT algorithm and the FT-PS mixed algorithm is still more than one order of magnitude faster than the PS algorithm alone. We also saw the GS algorithm helped improve the result of the PS algorithm; however, the GS algorithm was found in our studies to fail to converge in some cases, especially on noisy images. Therefore, we used the FT-PS mixed algorithm in the following studies.

One may wonder why the FT algorithm, based on a more precise formula, yields not as good results as the GS algorithm. The reason lies in the numerical discrete Fourier transform algorithm, which assumes periodical boundary conditions. Since in real experiments, the periodical boundary condition is generally not satisfied, the assumption may lead to deviation in the resultant phase maps. The GS algorithm, on the contrary, is based on Dirichlet, Neumann or mixed boundary conditions, thus may achieve better results than the FT algorithm.

Then we considered the robustness of the algorithms against image noise. We simulated the effect of image noise by adding a random number to each pixel value of the simulated images. Since the Poisson distributed shot noise can be approximated by a Gaussian distribution when the photon flux is high, we used a Gaussian distributed random number with a mean value zero and a standard deviation equal to a percentage of the pixel value. Thus a higher percentage means higher noise level.

In this test, we also tried the de-noise algorithm for removing Poisson noise presented by [21]. In this algorithm, the noise reduction is controlled by the weight of the fidelity term. The larger the weight is, the more details are kept; the smaller the weight is, the more smooth the resultant image is. In Figure 4(a), we show the SDDI results of the FT and the mixed algorithms on a set of images containing 1% noise and denoised using different fidelity weight. We found that in all cases the mixed algorithm gave better result than the FT algorithm for noisy images. The result also proved that the optimal range of the fidelity weight is from 0.1 and 0.3. Figure 4(b) shows the result of the mixed algorithm versus the noise level. One can see that the algorithm is quite robust for image noise below the 1% level. Evident deviation can be seen for noise level higher than 1%.

We tested the mixed algorithm experimentally on a dual detector x-ray phase imaging system [22]. The micro-focus tube (XTG Ultrabright Micro-focus Source, Oxford Instruments, Inc., Scotts Valley, CA) was operated at 40kVp. The average photon energy is measured as 19.52keV. The object-to-image distance *R*_{2} is 1.220m, resulting in a geometrical magnification of *M* = 1.66. Images were recorded with two mammographic CR detectors (Konica Minolta Medical Imaging USA, Inc., Wayne, NJ) with resolution of 43.75*μ*m. The sample object is a piece of chicken breast muscle scattered with several small pieces of calcium tablets. The muscle and calcium pieces resemble soft tissue and calcifications in mammography.

In the experiment, both the attenuation based image and the phase contrast image were obtained in one exposure using the dual detector scheme [23, 22]. We assumed that the thickness of the detector-1 used in the experiment is uniform, so the attenuation and phase maps of detector-1 are constants. Figure 5 shows (a) the attenuation based image and the phase retrieval result obtained by (b) the FT algorithm and (c) the mixed algorithm. One can see that both the FT algorithm and the mixed algorithm are capable of revealing the structure of the muscle and the calcium pieces. However, as can be seen from the background part of the images, the result by the FT algorithm is evidently hampered by the “periodical boundary condition” problem discussed above; while the mixed algorithm can in an extent improve the result.

We studied three iterative phase retrieval algorithms, based on Fourier transform, on Poisson solver and on the Gerchberg-Saxton algorithm, respectively. We found the FT algorithm is faster and the PS algorithm is more precise; the GS algorithm is much slower than the two and fails to converge in some cases. Therefore we suggest the FT-PS mixed algorithm for both higher efficiency and higher precision. The robustness of the mixed algorithm was tested using computer simulated noisy images. We tested the mixed algorithm experimentally by applying it to retrieve the phase map of a piece of chicken breast muscle.

The research is supported in part by the University of Oklahoma VPR office, and the Charles and Jean Smith Chair endowment fund.

**Publisher's Disclaimer: **This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Fanbo Meng, Center for Bioengineering and School of Electrical and Computer Engineering, University of Oklahoma, Norman, OK 73019, USA.

Da Zhang, Center for Bioengineering and School of Electrical and Computer Engineering, University of Oklahoma, Norman, OK 73019, USA.

Xizeng Wu, Department of Radiology, University of Alabama at Birmingham, Birmingham, AL 35233, USA.

Hong Liu, Center for Bioengineering and School of Electrical and Computer Engineering, University of Oklahoma, Norman, OK 73019, USA.

1. Wu X, Liu H. Clarification of aspects in in-line phase-sensitive x-ray imaging. Med Phys. 2007;34(2):737–743. [PubMed]

2. Meng FB, Liu H, Wu X. An iterative phase retrieval algorithm for in-line x-ray phase imaging. Opt Express. 2007;15(13):8383–8390. [PubMed]

3. Guigay JP, Langer M, Boistel R, Cloetens P. Mixed transfer function and transport of intensity approach for phase retrieval in the fresnel region. Opt Lett. 2007;32(12):1617–1619. [PubMed]

4. Gerchberg RW, Saxton WO. A practical algorithm for the determination of phase from image and diffraction plane pictures. Optik. 1972;35(2):237–246.

5. A Sommerfeld. Lectures on Theoretical Physics. Vol. 4. Optics, New York: Academic; 1954.

6. Lalor Éamon. Conditions for the validity of the angular spectrum of plane waves. J Opt Soc Am. 1968;58(9):1235.

7. Born M, Wolf E. Principles of Optics. 6th. Pergamon, Oxford: 1980.

8. Levi A, Stark H. Image restoration by the method of generalized projections with application to restoration from magnitude. J Opt Soc Am A. 1984;1:932–943.

9. Fienup JR. Phase retrieval algorithms: a comparison. Appl Opt. 1982;21(15):2758. [PubMed]

10. Marchesini S. A unified evaluation of iterative projection algorithms for phase retrieval. Rev Sci Instrum. 2007;78(1):011301. [PubMed]

11. Gureyev TE. Composite techniques for phase retrieval in the fresnel region. Opt Commun. 2003;220:49–58.

12. Zhang Y, Pedrini G, Osten W, Tiziani HJ. Phase retrieval microscopy for quantitative phase-contrast imaging. Optik. 2004;115(2):94–96.

13. Quiney HM, Nugent KA, Peele AG. Iterative image reconstruction algorithms using wave-front intensity and phase variation. Opt Lett. 2005;30(13):1638–1640. [PubMed]

14. Wu X, Liu H. A general theoretical formalism for x-ray phase contrast imaging. J X-ray Sci Tech. 2003;11(1):33–42. [PubMed]

15. Wu X, Liu H. Phase-space formulation for phase-contrast x-ray imaging. Appl Opt. 2005;44(28):5847–5854. [PubMed]

16. Gureyev TE, Nesterets YL, Paganin DM, Pogany A, Wilkins SW. Linear algorithms for phase retrieval in the Fresnel region. 2. Partially coherent illumination. Opt Commun. 2006;259(2):569–580.

17. Gureyev TE, Roberts A, Nugent KA. Partially coherent fields, the transport-of-intensity equation, and phase uniqueness. J Opt Soc Am A. 1995;12(9):1942.

18. Nugent KA, Gureyev TE, Cookson DF. Quantitative phase imaging using hard x rays. Phys Rev Lett. 1996;77(14):2961–2964. [PubMed]

19. Trescott PC, Larson SP. Comparison of iterative methods of solving twodimensional groundwater flow equations. wat. resourc. Wat Resour Res. 1977;13(1):125–136.

20. Douglas J. Alternating direction methods for three space variables. Numer Math. 1962;4:41–63.

21. Le T, Chartrand R, Asaki TJ. A variational approach to reconstructing images corrupted by poisson noise. J Math Imaging Vis. 2007;27(3):257–263.

22. Meng FB, Yan A, Zhou G, Wu X, Liu H. Development of a dual-detector x-ray imaging system for phase retrieval study. Nucl Instrum Meth B. 2007;254(2):300–306.

23. Wu X, Liu H. A dual detector approach for x-ray attenuation and phase imaging. J X-ray Sci Tech. 2004;12(1):35–42.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |