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

**|**HHS Author Manuscripts**|**PMC2925321

Formats

Article sections

- Abstract
- 1. Introduction
- 2. The attenuation-partition based algorithm (APBA) and an image accuracy measure
- 3. Simulation Tests
- 4. Discussion and Conclusion
- References and links

Authors

Related links

Opt Express. Author manuscript; available in PMC 2010 August 23.

Published in final edited form as:

PMCID: PMC2925321

NIHMSID: NIHMS226357

The publisher's final edited version of this article is available at Opt Express

See other articles in PMC that cite the published article.

The phase retrieval is an important task in x-ray phase contrast imaging. The robustness of phase retrieval is especially important for potential medical imaging applications such as phase contrast mammography. Recently the authors developed an iterative phase retrieval algorithm, the attenuation-partition based algorithm, for the phase retrieval in inline phase-contrast imaging [1]. Applied to experimental images, the algorithm was proven to be fast and robust. However, a quantitative analysis of the performance of this new algorithm is desirable. In this work, we systematically compared the performance of this algorithm with other two widely used phase retrieval algorithms, namely the Gerchberg-Saxton (GS) algorithm and the Transport of Intensity Equation (TIE) algorithm. The systematical comparison is conducted by analyzing phase retrieval performances with a digital breast specimen model. We show that the proposed algorithm converges faster than the GS algorithm in the Fresnel diffraction regime, and is more robust against image noise than the TIE algorithm. These results suggest the significance of the proposed algorithm for future medical applications with the x-ray phase contrast imaging technique.

Differing from the conventional x-ray imaging, which relies on the small differences in x-ray attenuation changes between tissues variable structure, inline phase contrast imaging is based on the tissues’ phase-shifts diffraction from the object to the detector. Since x-ray phase-shift differences between tissue and lesions are about one thousand times larger than attenuation differences [2, 3, 4], x-ray phase contrast imaging has the potential to enhance the lesion detection sensitivity and specificity, and reduce the radiation dose compared to conventional x-ray imaging. In the inline phase contrast imaging, the diffracted phase-shifts form bright and dark fringes at tissue boundaries and this bright and dark fringes are called edge enhancement. The edge enhancement relies on the spatial coherence of the x-ray source, the Laplacian and gradients of x-ray phase-shifts caused by the tissue, and the gradients of the objects attenuation [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]. One procedure of phase contrast imaging is to disentangle tissue phase-shifts from the mixed contrast mechanism and recover the phase maps from acquired phase contrast images. This procedure is called phase retrieval. Phase retrieval technique plays a central role in phase contrast x-ray imaging. By means of phase retrieval, one can reconstruct a quantitative map of phase-shifts, a phase image of the imaged object [4, 7, 14, 16, 17]. The amount of x-ray phase-shifts by tissues is determined by

$$\phi (\overrightarrow{\mathbf{r}})=-\left(\frac{hc}{E}\right){r}_{\text{e}}\int {\rho}_{e}(\overrightarrow{\mathbf{r}},z)dz=-\left(\frac{hc}{E}\right){r}_{\text{e}}{\rho}_{\text{e},\text{p}}(\overrightarrow{\mathbf{r}}),$$

(1)

where *r*_{e} is the classical electron radius, *h* the Plank constant, *c* the speed of light, *E* the x-ray photon energy, and *ρ*_{e,p}, the integration of the electron density *ρ _{e}* over the x-ray path, is called the projected electron density [2, 3, 4]. So a retrieved phase map is equivalently a map of imaged object’s quantitative projected electron densities. Moreover, phase retrieval is also a necessary procedure for phase-sensitive volumetric imaging, such as the phase sensitive tomography and tomosynthesis, to acquire the artifact free 3D images of object attenuation coefficients and electron densities [8, 15, 16].

Phase retrieval is based on x-ray propagation equation derived either from the Fresnel diffraction or the Wigner distribution based phase-space formalism [5, 18, 9, 19, 20]. To be specific, let (**$\stackrel{\u20d7}{r}$**) be the x-ray phase-shift caused by the imaged object, and *A*_{0}(**$\stackrel{\u20d7}{r}$**) the x-ray transmission, or the attenuation-map of the object. Then the Fourier transform of the x-ray intensity, (*I*)(*$\stackrel{\u20d7}{u}$*) = *∫*_{2} *I*(*$\stackrel{\u20d7}{x}$*)exp [2 *πi$\stackrel{\u20d7}{x}$* · *$\stackrel{\u20d7}{u}$*]*d$\stackrel{\u20d7}{x}$*, at position away from the object with distance *R*_{2}, of the monochromatic point x-ray source starting at a place away from the object *R*_{1} distance, can be modeled by the following [9]

(2)

where *I*_{in} is the incident x-ray intensity at *R*_{1}, *λ* the wavelength of the monochromatic point x-ray source and *M* = (*R*_{1} + *R*_{2})/*R*_{1} is the geometric magnification. When the FT-Space Fresnel propagator *πλ R*_{2}*$\stackrel{\u20d7}{u}$*^{2}/*M* 1, Eq. (2) can be simplified to the Transport of Intensity Equation (TIE) [21, 4, 9]

$$I(\overrightarrow{\mathbf{r}};{R}_{1}+{R}_{2})=\frac{{I}_{\text{in}}}{{M}^{2}}\left\{{A}_{0}^{2}\left(\frac{\overrightarrow{\mathbf{r}}}{M}\right)-\frac{\lambda {R}_{2}}{2\pi M}\nabla \xb7({A}_{0}^{2}\nabla \phi )\left(\frac{\overrightarrow{\mathbf{r}}}{M}\right)\right\}.$$

(3)

It is worthy to note that Eq. (3) is valid only for moderate resolution images. For high resolution images, i.e. when the FT-Space Fresnel propagator *πλ R*_{2}*$\stackrel{\u20d7}{u}$*^{2}/*M* is close or greater than *π*, any phase retrieval algorithms based on Eq. (3) need to be substituted to Eq. (2) [22, 23, 24]. In this paper, the algorithms discussed are all based on Eq. (3), i.e. the moderate resolution is satisfied.

Numerous algorithms have been suggested on how to effectively recover the phase-shift from the phase contrast images. Among these, two algorithms are most widely used. One is the TIE algorithm implemented by Allen and Oxley in [25]. The other is the GS algorithm developed first by Gerchberg and Saxton in [26] and later improved by Fienup [27, 28]. These two algorithms have both their advantages and disadvantages. The TIE algorithm is a direct approximate method which is fast but is unstable with noisy data; the GS algorithm on the other hand is relatively more stable than the TIE algorithm[25] but the converging speed is slow especially for the field of medical imaging. In [1], the authors developed an Attenuation-Partition Based Algorithm (APBA) based on the phase-attenuation duality property of soft tissues under higher x-ray energy. This algorithm is fast and stable for potential medical imaging. We compared the performance of this algorithm with the TIE algorithm for two groups of data under the condition of medical imaging in [1], including the phase retrieval from phase-contrast images of a breast lumpectomy specimen. In this paper, we will make a systematic analysis about this algorithm and compare its performance with the one of the GS algorithm and the TIE algorithm with simulated data.

The paper is organized as follows. In Section 2, we first summarize the attenuation-partition based algorithm (APBA), which is motivated by our observation of the phase-attenuation duality[14]. Then we give a measure, called total variation, used to evaluate the closeness of two image data. This measure is used as a quantitative standard in comparing the performance between different algorithms in the following section. In Section 3, we first develop a breast specimen model which can reflect the attenuation and phase changes with respect to the x-ray energy change (Section 3.1), and then compare the performance of the algorithm with the GS algorithm (Section 3.2) and the TIE algorithm (Section 3.3). Finally, we conclude this paper in Section 4.

The attenuation-partition based algorithm (APBA) is a recently developed iterative algorithm for phase retrieval[1]. It was derived from our previous notion of the phase-attenuation duality[14], and it takes advantage of the correlation between the attenuation and phase-shift for phase retrieval. As is well known, tissue’s attenuation change *A*_{0} in the diagnostic x-ray imaging arose from three x-ray and tissue interactions: the photoelectric absorption, the coherent scattering, and the incoherent scattering[29, 30, 9, 14]. However, among the three interactions, the attenuation caused by incoherent scattering *A*_{KN}, which is dominated by the x-ray Compton scattering, deserves a special attention. This is because both of *A*_{KN} and the x-ray phase-shift are determined by the tissue’s projected electron density:

$${A}_{\text{KN}}(\overrightarrow{\mathbf{r}})=exp\left[-\frac{{\sigma}_{\text{KN}}}{2}{\rho}_{\text{e},\text{p}}(\overrightarrow{\mathbf{r}})\right],\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phi (\overrightarrow{\mathbf{r}})=-\lambda {r}_{\text{e}}{\rho}_{\text{e},\text{p}},$$

(4)

where *λ* is the x-ray wavelength, *r*_{e} the classical electron radius, *ρ*_{e,p} the projected electron density as defined in Eq. (1), and *σ*_{KN} is the total cross section for x-ray Compton scattering with a free electron:

$${\sigma}_{\text{KN}}({E}_{\text{photon}})=2\pi {r}_{\text{e}}^{2}\left\{\frac{1+\eta}{{\eta}^{2}}\left[\frac{2(1+\eta )}{1+2\eta}-\frac{1}{\eta}log(1+2\eta )\right]+\frac{1}{2\eta}log(1+2\eta )-\frac{1+3\eta}{{(1+2\eta )}^{2}}\right\},$$

(5)

with *η* = *E*_{photon}/*m*_{e}*c*^{2}. Here we denote the photon energy of the primary x-ray beam by *E*_{photon} and *m*_{e}*c*^{2} is the rest electron energy. Eq.(4) suggests clearly that the x-ray attenuation and phase-shift by tissue has certain correlation. Our idea is to utilize this correlation for facilitating the phase retrieval.

Of course, the extent of this correlation between phase and attenuation depends on the x-rays photon energy as well as the tissues physical composition. For example, for light elements with atomic numbers *Z* ≤ 10, x-ray attenuation is dominated by the Compton scattering for x-rays of 60 keV or higher, i.e. *A*_{0} ≈ *A*_{KN}[14]. We call this relationship between phase-shift and attenuation the phase-attenuation duality. The phase-attenuation duality can be used for phase retrieval as follows. Consider a phase contrast imaging setting with a point source of wavelength *λ*. The object is at a distance *R*_{1} from the source. We denote *R*_{2} as the distance from object to detector plane, *M* = (*R*_{1} + *R*_{2})/*R*_{1} the geometric magnification factor, *I*_{in} the entrance x-ray intensity at *R*_{1}, and *I _{D}*(

(6)

where

$$\stackrel{\sim}{k}=\frac{\lambda {R}_{2}}{2\pi M}\xb7\frac{\lambda {r}_{\text{e}}}{{\sigma}_{\text{KN}}},$$

(7)

and , for sake of convenience, is called the “duality transform” acting on the normalized image *I*.

In general imaging cases, such as with low energy x-rays or an object contains calcified tissues such as calcification, this phase-attenuation duality does not hold. However, we can still factor out tissue’s total attenuation *A*_{0} as

$${A}_{0}(\overrightarrow{\mathbf{r}})={A}_{\text{KN}}(\overrightarrow{\mathbf{r}})\xb7{A}_{\text{pe},\text{coh}}(\overrightarrow{\mathbf{r}}),$$

(8)

where we denote the attenuation caused by photoelectric absorption and coherent scattering by *A*_{pe,coh}. Strictly speaking, *σ*_{KN} is only Compton scattering cross-section, it may be slightly different from the incoherent scattering cross-section for high-Z elements. This is because while Compton scattering is x-ray scattering from the free electrons, the incoherent scattering is that from the bound atomic electrons[31]. So when we factor *A*_{0} = *A*_{KN} · *A*_{pe,coh} (Eq. (8)) we actually factor the small residual binding effect of atomic electrons into *A*_{pe,coh}. With this understanding, Eqs. (4) and (8) are rigorously valid. The notion of Eq. (6) and Eq. (8) led us to the development of the attenuation-partition based algorithm [1]. While the derivations and the algorithm details of this method can be found from Ref. [1], a brief description of the method is as follows. Our goal is to retrieve the phase map from the two normalized images: one is the object’s attenuation image
${A}_{0}^{2}$ acquired with *R*_{2} = 0, and the other is the acquired phase contrast image *I* = *M*^{2} · *I _{D}*(

$${A}_{0}={A}_{\text{KN}}-\delta A,\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\delta A={A}_{\text{KN}}(1-{A}_{\text{pe},\text{coh}}),$$

(9)

and find the correction term *δ A* using the estimate of *A*_{KN}. We then employing the Fresnel Diffraction transform (as defined in Eq. (11)) to transport the wave field *δ Ae ^{i}* from

The above iterative procedure can be summarized in flow chart Fig. 1 and the Algorithm

In an imaging experiment, assume we have performed two imaging measurements. One image (radiograph) is the attenuation image
${A}_{0}^{2}$ acquired at SID= R_{1}, another is a normalized phase-contrast image I = M^{2} · I_{D}($\stackrel{\u20d7}{r}$_{D})/I_{in} acquired at SID= R_{1} + R_{2}. With
${A}_{0}^{2}$ and I as well as the initial δ I, usually 0, as known inputs, we first assume
${I}_{\text{KN}}={\left(\sqrt{I}+\sqrt{\delta I}\right)}^{2}$. Then

- Compute ${A}_{\text{KN}}=\sqrt{\mathfrak{D}({I}_{\text{KN}})}$ and from the duality transform Eq. (6).
- Compute δ A from$$\delta A={A}_{\text{KN}}\xb7(1-P),\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}P={A}_{0}/{A}_{\text{KN}}.$$(10)Equations (9) and (10) are in fact the same equations. The advantage of Eq. (10) over (9) is that we can set a threshold for P. We know P = A
_{pe,coh}in the ideal case and A_{pe,coh}is bounded between 1 and A_{0}. The computation rounding error or the presence of measured noise in the acquired data can make the value of P over pass these bounds in the iterative computations. By setting a threshold upper bound ubd = 1 and lower bound lbd = min(A_{0}), the minimum value of A_{0}, to P in the iterative computations, we can make the algorithm more stable. Moreover, if we know a better lbd for A_{pe,coh}, other than the minimum of A_{0}, the converging speed of the algorithm can be greatly improved. - Compute δI by Fresnel propagate δAe
^{i}from position R_{1}to R_{2}: δI = | (δAe^{i}|^{2}with$${\mathfrak{F}}_{\mathfrak{r}}(T)(\overrightarrow{r})=\frac{1}{\lambda {R}_{2}}{\int}_{{\mathbb{R}}^{2}}exp\left[i\frac{\pi M}{\lambda {R}_{2}}{\left(\frac{\overrightarrow{r}}{M}-\overrightarrow{\xi}\right)}^{2}\right]T(\overrightarrow{\xi})d\overrightarrow{\xi}.$$(11) - Compute ${I}_{\text{KN}}={\left(\sqrt{I}+\sqrt{\delta I}\right)}^{2}$. Go to (1) for next iteration.

The number of iterations or an accuracy measure can be used to determine when to exit the program: assuming ||·|| is some kind of norm, that can effectively reflect the accuracy of the retrieved data as an image, if ||*δI _{k}*

An appropriate image accuracy measure should be a measure that can effectively reflect the accuracy of the data AND at the same time correctly reflect the visual perception of the data as an image, since an image’s visual perception is crucial for diagnostic radiology. In the next subsection the authors suggest a measure which can be employed as an appropriate image accuracy measure, which was first investigated by Rudin in [32].

A continuous signal is generally represented as a function of vector variables: *f*(**$\stackrel{\u20d7}{r}$**). A sampled signal will be represented as a one- (or higher) dimensional sequence of real numbers. In this paper, we will denote the continuous two dimensional image as an intensity function of two dimensional variables, such as *f* (*x*,*y*), or *f* (**$\stackrel{\u20d7}{r}$**), **$\stackrel{\u20d7}{r}$** = (*x*,*y*). The sampled 2-D image will be represented by *f* (*i*, *j*), *i* = 1,2, ···, *m*, *j* = 1,2, ···, *n*. In practice, to estimate a true signal in noise, the most frequently used methods are based on the least squares criteria and thus an intuitive measure for the closeness of two image functions *f* and *g* is, similar to statistical standard deviation, based on:

$$\text{std}(g,f):=\frac{{\left[{\int}_{\mathrm{\Omega}}{(g(\overrightarrow{\mathbf{r}})-f(\overrightarrow{\mathbf{r}})-\mu )}^{2}d\overrightarrow{\mathbf{r}}\right]}^{1/2}}{V(\mathrm{\Omega})},$$

(12)

where Ω is the finite domain of image functions *f* and *g*, *V* (Ω) represents the area of domain Ω, and *μ* = *∫*_{Ω}(*g* − *f)d***$\stackrel{\u20d7}{r}$**/*V* (Ω) is the statistical mean value of *g* − *f*. In [32], L. I. Rudin investigated the relation of edge formation of the 2-D digital image and the singularities of the image function and pointed out that the image intensity function belongs to the space of functions of bounded total variation. Rudin et al.[33] pointed out that the proper norm for images is the total variation (TV) norm, which is the *L*_{1} norm of the gradient of the image function, and not the *L*_{2} norm. For two image functions *f* and *g*, we define the TV norm of *g* − *f* as the closeness of the two images:

$$\text{TV}(g,f):=\frac{{\int}_{\mathrm{\Omega}}\mid \nabla (g-f)\mid d\overrightarrow{\mathbf{r}}}{V(\mathrm{\Omega})}=\frac{{\int}_{\mathrm{\Omega}}\sqrt{{\left({\scriptstyle \frac{\partial (g-f)}{\partial x}}\right)}^{2}+{\left({\scriptstyle \frac{\partial (g-f)}{\partial y}}\right)}^{2}}d\overrightarrow{\mathbf{r}}}{V(\mathrm{\Omega})},$$

(13)

where is the gradient operator. Since std is the form of *L*_{2} norm, it is not suitable as a measure to represent the closeness between two images. For example, for the three “Lena” image functions shown in Fig. 2, the std measures of _{1} − _{2} and _{1} − _{3} are std(_{1}, _{2}) = 0.206, std(_{1}, _{3}) = 0.472 respectively. But the image _{3} is much more closer to _{1} than _{2} does in visual perception. This is because the digital representation of an image depends not only on the pixel values, it also depends on the contrast changes between neighbor pixels. An image’s visual perception is crucial for diagnostic radiology. This contrast changes between neighbor pixels can better be represented by the gradient changes of the image functions. For example, the TV measures of _{1} − _{2} and _{1} − _{3} are TV(_{1}, _{2}) = 0.0247 and TV(_{1}, _{3}) = 0.0138 respectively, more appropriate in reflecting the visual perception. In this paper, we will use the TV norm (13) as the measure for closeness between two compared image functions.

Note that the TV norm between two image functions, say *g* and *f*, equals 0 if and only if *g* differs from *f* by a constant. This feature does not affect its appropriateness for phase retrieval since it is well known that the recovered phase is unique up to a constant with given information about the attenuation-map *A*_{0} and the phase-contrast intensity-map *I*.

In order to investigate the performance of the algorithm constructed above, we perform computer simulations in this section. In Section 3.1, we first construct a breast tissue model that represents the phase-shifts and attenuation of breast tissues and embedded microcalcifications for different x-ray energies. In our simulation tests, the distances of source point to object, *R*_{1}, and object to detector, *R*_{2}, are set to 26 inches (0.66 m) respectively. In this way the magnification factor *M* = (*R*_{1} + *R*_{2})/*R*_{1} is equal to 2. For convenience, the incident x-ray intensity *I*_{in} at *R*_{1} is set to *M*^{2} (4). We compare the performance of our algorithm, APBA, with two other widely used algorithms: the Gerchberg-Saxton (GS) algorithm (Section 3.2) and the TIE algorithm (Section 3.3).

The tissue’s phase-shift and attenuation are determined not only by the tissue’s physical composition, but also by the x-ray energy. With different x-ray energy, the same tissue has different phase-shift and attenuation change. Simulation models used in literature often do not incorporate these changes. In this subsection, we construct a breast specimen model that can represent tissue attenuation and phase shifts according to employed x-ray energies as well as tissue’s compositions.

In our model the tissue has two physical compositions: the 50% glandular-50% adipose breast tissue and the microcalcifications. In order to simulate the morphological aspects of breast tissues, we adopted a radiograph of a breast lumpectomy specimen with pixel values rescaled and the metal localization wire removed by replacing the pixel value at wire position with a mean pixel value at near by positions. It is difficult to remove all the residual trace of the wire this way as can be seen from the following image display especially for the phase contrast image (Fig. 5(c)). In the phase contrast image (Fig. 5(c)), the small residual variation from the original wire-track really got enhanced. The linear attenuation coefficients and electron densities for the 50% glandular-50% adipose breast tissues are computed from the tissue’s elemental composition and the interpolated elemental mass attenuation coefficients from the standard reference in the mammographic radiation dosimetry [34, 35]. Moreover, each mass attenuation coefficient is further broken down to two components: one from x-ray photoelectric absorption and coherent scattering, and another from incoherent scattering. In this way, the total attenuation is partitioned as a product of ${A}_{\text{pe},\text{coh}}^{2}$ and ${A}_{\text{KN}}^{2}$ as defined in Eq. (8) above. In addition, to simulate the microcalcifications in breast, four small ellipsoids of calcium carbonate (CaCO3) are embedded in the specimen model. The diameter of the ellipsoids can be adjusted in simulating different size of the microcalcifications. In the following simulations, to test phase-contrast sensitivity, we set the diameters of the four ellipsoids to 10, 5, 10, and 5 microns in x-ray direction and 300, 200, 300, 200 microns in detector plane respectively.

Image representation of the inputs generated from the simulation model and Fresnel propagation. (a) the phase map ; (b) the attenuation map
${A}_{0}^{2}$; and (c) the normalized Fresnel propagated phase contrast image *I* with object to detector distance **...**

The attenuation image ${A}_{0}^{2}$ of the specimen model simulated with 35.5 keV x-ray, and its two corresponding partition images ${A}_{\text{KN}}^{2}$ and ${A}_{\text{pe},\text{coh}}^{2}$, are shown in Fig. 3. For a comparison, the profiles of ${A}_{\text{KN}}^{2}$ and ${A}_{\text{pe},\text{coh}}^{2}$ along the line passing through the microcalcifications are shown in Fig. 4 simulated with x-ray energy equals 18.5, 35.5 and 59.5 keV, respectively. We can see that the contribution of ${A}_{\text{pe},\text{coh}}^{2}$ to the total attenuation ${A}_{0}^{2}$ gets smaller when x-ray energy is getting higher. Especially, when x-ray energy equals 59.5 keV, the contribution of ${A}_{\text{pe},\text{coh}}^{2}$ for the soft tissue can almost be neglect, as is expected.

Image manifest of (a)
${A}_{\text{pe},\text{coh}}^{2}$, (b)
${A}_{\text{KN}}^{2}$ and (c)
${A}_{\text{0}}^{2}$ when x-ray energy equals 35.5 keV.

The GS algorithm is an iterative algorithm for phase retrievals from a pair of images at two planes related by the Fourier transform. For details readers are referred to [26]. The GS algorithm is a classical algorithm which is widely used in electron microscopy, wave front sensing, astronomy, crystallography, and many other fields involving phase recovery [27, 28, 36].

By replacing the Fourier Transform in the GS algorithm with the Fresnel transform of Eq. (11), we get a modified version of the GS algorithm extended to the *Fresnel diffraction regime*. Our previous simulation tests showed that this modified GS algorithm converges very slow for object-detector distance *R*_{2} ≈ 1 m, and converges faster for larger *R*_{2}, such as images acquired at synchrotron beam lines. But this is generally not applicable to the field of clinical imaging, due to the physical constraints such as compact sizes of hospital rooms.

In our simulation tests, we compare the performance of APBA and that of the GS algorithm. The photon energy of the point x-ray source is set to 35.5 keV, and the distances from the source to object and object to detector are set to *R*_{1} = *R*_{2} = 26 inches (0.66 m) respectively. For convenience, the incident x-ray intensity *I*_{in} at *R*_{1} is set to *M*^{2}, where *M* = (*R*_{1} + *R*_{2})/*R*_{1} is the magnification factor. The phase map and the attenuation
${A}_{0}^{2}$ are generated from our phantom model for 35.5 keV x-ray. Fig. 5 shows the simulated images of the phase map , the attenuation image
${A}_{0}^{2}$ and the phase-contrast image *I*.

In the simulation test, the iteration for our attenuation-partition based algorithm (APBA) and the GS algorithm are performed 100 steps. The corresponding recovered phase images are shown in Fig. 6(b) and (c). In Fig. 6(a), the solid line represents the change of total variation (TV) of the retrieved phase using attenuation-partition based algorithm with respect to the iteration steps. The dashed line represents the change of TV of the retrieved phase using the (modified) GS algorithm with respect to iteration steps. We can see that the converging speed of the (modified) GS algorithm is much slower than that of attenuation-partition based algorithm (the change of TV of APBA from step 1 to step 100 is 1.04*E* – 3, almost 33 times greater than that, 3.17*E* – 5, of the GS algorithm.) In addition, from the visual perception point of view, we see that the phase map retrieved with the attenuation-partition based algorithm (APBA) is much better than that retrieved with the (modified) GS algorithm.

Comparison of the performance of the GS algorithm and APBA. (a) plot of the accuracy measures with respect to iteration steps. The plot with solid line represents the APBA. The one with dashed line represents the GS algorithm; (b) recovered phase map **...**

The main difference between APBA and the GS algorithm is that APBA takes the advantage of the phase-attenuation correlation property, although the extent of this correlation is not known in *priori*, but the GS algorithm does not. From this example we see that the phase-attenuation correlation is a very important information that shouldn’t be neglected in the algorithm development for phase retrieval.

We have mentioned the transport of intensity equation in Eq. (3) in Section 1. Since Teague proposed the TIE algorithm for phase retrieval in 1983 [21], numerous phase retrieval algorithms have been suggested on how to effectively search the numerical solution of the TIE [4, 20, 25, 37, 38, 39, 40]. Among the methods of solving the above TIE for the phase map, the one based on Fast Fourier Transform, proposed by Nugent et al. [4], and Allen and Oxley in [25], is the most widely used. In the form of pseudo-differential operator, the solution phase map is given by

$$\phi (\overrightarrow{\mathbf{r}})=-\frac{2\pi M}{\lambda {R}_{2}}{\nabla}^{-2}\left\{\nabla \xb7\left[\frac{\nabla [{\nabla}^{-2}(I-{A}_{0}^{2})]}{{A}_{0}^{2}}\right]\right\},$$

(14)

for the given normalized phase-contrast image *I* and attenuation image
${A}_{0}^{2}$. Here is the gradient operator, and ^{−2} is the inverse Laplacian operator.

The advantage of this algorithm is that it does not require the boundary information in solving the TIE (assuming the image data is periodic); it is a deterministic method and thus the algorithm is fast and accurate comparing to most iterative algorithms. In this section we compare the performance of the TIE algorithm and that of APBA for two kind of cases: first for the ideal case without any noise and any image misalignment, then for cases simulating practical situation with x-ray imaging noise and possible image misalignment. In these simulation tests, the imaging geometries are the same as in the previous subsection, and x-ray energy is again 35.5 keV.

For the ideal case without any noise and any image misalignment, the performance comparison results are shown in Fig. 7. For the ideal case the TIE algorithm is accurate both in TV measure and visual perception. APBA can also achieves this accuracy but needs 1110 iteration steps to have its TV measure of 0.00215226, an error smaller than that of 0.00215269 with the TIE algorithm.

Comparison of APBA and the TIE algorithm with pure data. (a) plot of the accuracy measure with respect to iteration steps. The plot with solid line represents the APBA. The one with dashed line represents the TIE algorithm. It needs 1110 steps for the **...**

However, the real test lies in the performance for the cases simulating practical situation with x-ray imaging noise and possible image misalignment. Obviously only the performance in these cases really matters in phase contrast imaging applications, especially for clinical imaging applications where the imposed radiation limits dictates existence of substantial x-ray quantum noise in acquired images. Implemented in the Fourier space, the inverse Laplacian ^{−2} in Eq. (14) is singular at zero spatial frequency. This singularity will amplify the noise in the images and result in failure of accurate phase retrieval for the TIE algorithm. To overcome this difficulty, some kind of regularization must be used. The most widely used regularization is Tikhonov regularization, which replaces ^{−2} by |**$\stackrel{\u20d7}{u}$**|^{2}/(|**$\stackrel{\u20d7}{u}$**|^{4} + *κ*^{2}), for some “favorable” constant parameter *κ*, called the Tikhonov regularization parameter. In this regularization scheme, the singularity is regularized, and the favorable parameter *κ* means the retrieved phase * _{κ}* corresponding this

$${x}_{\kappa}=arg\underset{x\in X}{min}{\left|\right|Ax-y\left|\right|}_{Y}^{2}+\kappa {\left|\right|x\left|\right|}_{X}^{2}.$$

(15)

It is *L*_{2} norm dependent. Since the proper norm for image data is the total variation (TV) norm [33], a favorable solution under the Tikhonov regularization principle can not be guaranteed a best solution in visual perception. Moreover, for relatively noisy acquired images, the TIE-based algorithm, even with Tikhonov regularization, often failed to retrieve the phase maps [1]. In the following, we will compare the performance of APBA and the TIE-algorithm when noise is present. In the practice of phase retrieval, there are generally two kinds of image data errors. One is the noise associated with image acquisitions, including the quantum noise of x-ray photon detections and detector electron noise. We assume the quantum noise dominates as is the case in most imaging applications. The other is the error caused by the misalignment between the attenuation map
${A}_{0}^{2}$ and phase-contrast image *I*. This is because usually the attenuation image and corresponding phase contrast image are generally acquired in two separate x-ray exposures. The misalignment could be resulted from the shift or tilting of the object or detector between the exposures. In the simulation, we associate each pixel value of an image a photon count *N*, so that *P*(*i*, *j*) = *c* · *N*(*i*, *j*), where *c* is a constant. Assuming the noise has a Poisson distribution, with variance *σ*^{2} = *N* at each pixel. In the simulation we assign a background noise level for each simulated image. This background noise level is defined as the ratio *δ _{b}* =

Comparison of the TIE algorithm and APBA with noise added. (a) True phase map ; (b) attenuation map
${A}_{0}^{2}$; (c) the normalized phase contrast image *I*; (d) recovered phase map with the TIE algorithm, no Tikhonov regularization is used; (e) recovered **...**

Three simulation tests are performed in the comparison: case 1: assume the acquired data
${A}_{0}^{2}$ and *I* has noise present but perfectly aligned; case 2: assume acquired data has no noise but has one pixel misalignment horizontally; case 3: assume acquired data has combined detector noise as well as one pixel misalignment. One bias in simulation for the TIE-algorithm should be mentioned: with the known true phase value, a favorable Tikhonov regularization parameter *κ* can be searched, but in practice this search is hardly feasible. In each case mentioned above, three phase retrievals are performed: 1. using the TIE-algorithm without Tikhonov regularization; 2. using the TIE-algorithm and the favorable Tikhonov regularization parameter; 3. using APBA with 10 step iteration. The Total variation (TV) of the results are listed in Tab. 1 and the recovered phase images for case 3 are shown in Fig. 8(d) – (f). The results using Tikhonov regularization are better than those not using Tikhonov regularization but worse than those using APBA. The influence of misalignment to APBA is little but disaster to the TIE algorithm. From the profile, shown in Fig. 9, along a line passing through the microcalcifications we can see that the values of the phase recovered from APBA is very close to the true phase value, but the values of the phase recovered from the TIE algorithm are distorted.

Profiles, along a line passing through the microcalcifications, of the recovered phase using APBA, the solid line, and using the TIE algorithm with Tikhonov regularization, the dashed line, in Case 3. The dash-dotted line is the true phase.

In above analysis, we assumed that the x-ray source is a quasi-monochromatic point source. In practice, one often employs conventional incoherent and polychromatic sources such as x-ray tubes for imaging. In the experiments with an x-ray tube source, the previous formula Eqs. (2–3) of the phase contrast image intensities should be modified. Since in the APBA method the duality transform is derived based on Eq. (3), hence the Eq. (6) should be modified accordingly. In our previous works we have studied this problem [20]. With the Wigner function based phase space formalism, we have proved that the coherence degree of a finite-size focal spot can be accounted for by the optical transfer function *OTF _{G.U}*.(

$${\mathit{OTF}}_{G.U.}\left(\frac{\overrightarrow{u}}{M}\right)=\frac{\int {I}_{\text{spot}}(\overrightarrow{\xi})exp\left[i2\pi \overrightarrow{\xi}\xb7{\scriptstyle \frac{(M-1)\overrightarrow{\mathbf{u}}}{M}}\right]d\overrightarrow{\xi}}{\int {I}_{\text{spot}}(\overrightarrow{\xi})d(\overrightarrow{\xi})}$$

where *I*_{spot}() is the intensity distribution of the focal spot. We found the generalized TIE equation with an x-ray tube source, that is, the x-ray intensity at the detector is given by [20]:

where operator denotes the convolution, *A*^{2} is the total attenuation of the imaged object, is the spectrally averaged phase-shift caused by the object, and · means the spectral average. Compare above equation with Eq. (3) and it is clear that the TIE-based phase retrieval method needs only two modifications: (i). Fourier space de-convolution of the measured intensity from *OTF _{G.U.}* (

$$\langle \stackrel{\sim}{k}\rangle =\frac{{r}_{e}{R}_{2}}{2\pi M}\xb7\langle \frac{{\lambda}^{2}}{{\sigma}_{\text{KN}}}\rangle .$$

Otherwise, the APBA flow chart is the same as that for the case with a quasi-monochromatic point source.

Phase retrieval is a crucial step for quantitative imaging such as reconstructing the 3-D distribution of tissue linear attenuation coefficients and refraction indices. However, images acquired in medical applications are relatively noisy, due to the radiation dose constraints, with low phase contrast effect, due to physical constraints such as compact sizes of hospital rooms. An phase retrieval algorithm which is robust to noise is necessary for potential medical phase contrast imaging. In [1], the authors developed an algorithm, called attenuation-partition based phase retrieval algorithm. It is an iterative algorithm which takes advantage of the correlation between the attenuation and phase-shift. The phase retrieval results from experimental images show that this algorithm is fast and robust [1]. In this work, we systematically compared the performance of this algorithm with other two widely used phase retrieval algorithms, namely the Gerchberg-Saxton (GS) algorithm and the Transport of Intensity Equation (TIE) algorithm. The systematical comparison is conducted by analyzing phase retrieval performances with a digital breast specimen model. We show that the proposed algorithm converges faster than the GS algorithm in the Fresnel diffraction regime, and is more robust against image noise than the TIE algorithm. These results suggest the significance of the proposed algorithm for future medical applications with the x-ray phase contrast imaging technique.

This research was supported in part by the Department of Defense Breast Cancer Research Program under award number W81XWH-08-1-0613 and the NIH grant 1R01CA142587.

1. Yan A, Wu X, Liu H. An attenuation-partition based iterative phase retrieval algorithm for in-line phase-contrast imaging. Opt Express. 2008;16:13,330–13,341. [PubMed]

2. Wilkins S, Gureyev T, Gao D, Pogany A, Stevenson A. Phase-contrast imaging using polychromatic hard X-rays. Nature. 1996;384:335–338.

3. Snigirev A, Snigireva I, Kohn V, Kuznetsov S, Shelokov I. On the possibilities of x-ray phase contrast microimaging by coherent high-energy synchrotron radiation. Rev Sci Instrum. 1995;66:5486–5492.

4. Nugent K, Gureyev T, Cookson D, Paganin D, Barnea Z. Quantitative Phase Imaging Using Hard X Rays. Phy Rev Lett. 1996;77:2961–2965. [PubMed]

5. Pogany A, Gao D, Wilkins S. Contrast and resolution in imaging with a microfocus x-ray source. Rev Sci Instrum. 1997;68:2774–2782.

6. Arfelli F, Bonvicini V, et al. Mammography with synchrotron radiation: phase-detected Techniques. Radiology. 2000;215:286–293. [PubMed]

7. Paganin D, Mayo S, Gureyev T, Miller P, Wilkins S. Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object. J Microsc. 2002;206:33–40. [PubMed]

8. Mayo S, Davis T, Gureyev T, Miller P, Poganin D, Pogany A, Stevenson A, Wilkins S. X-ray phase-contrast microscopy and microtomography. Opt Express. 2003;11:2289–2302. [PubMed]

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

10. Wu X, Liu H. Clinical implementation of phase-contrast x-ray imaging: Theoretical foundations and design considerations. Med Phys. 2003;30:2169–2179. [PubMed]

11. Wu X, Liu H. A new theory of phase-contrast x-ray imaging based on Wigner distributions. Med Phys. 2004;31:2378–2384. [PubMed]

12. Donnelly E, Price R, Pickens D. Experimental validation of the Wigner distributions theory of phase-contrast imaging. Med Phys. 2005;32:928–931. [PubMed]

13. Zhang D, Donvan M, Fajardo L, Archer A, Wu X, Liu H. Preliminary feasibility study of an in-line phase contrast x-ray imaging prototype. IEEE Trans Biomed Eng. 2008;55:2249–2257. [PubMed]

14. Wu X, Liu H, Yan A. X-ray phase-attenuation duality and phase retrieval. Opt Lett. 2005;30(4):379–381. [PubMed]

15. Wu X, Liu H. X-Ray cone-beam phase tomography formulas based on phase-attenuation duality. Opt Express. 2005;13:6000–6014. [PubMed]

16. Cloetens P, Mache R, Schlenker M, Lerbs-Mache S. Quantitative phase tomography of Arabidopsis seeds reveals intercellular void network. PNAS. 2006;103:14,626–14,630. [PubMed]

17. Wu X, Liu H, Yan A. Phase-Contrast X-Ray Tomography: Contrast Mechanism and Roles of Phase Retrieval. Eur J Radiology. 2008;68:S8–S12. [PubMed]

18. Paganin D, Nugent K. Noninterferometric Phase Imaging with Partially Coherent Light. Phy Rev Lett. 1998;80:2586–2589.

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

20. Wu X, Liu H. Phase-space evolution of x-ray coherence in phase-sensitive imaging. Appl Opt. 2008;47:E44–E52. [PubMed]

21. Teague M. Deterministic phase retrieval: a Green’s function solution. J Opt Soc Am. 1983;73:1434–1441.

22. Gureyev T, Nesterets Y, Paganin D, Pogany A, Wilkins S. Linear algorithms for phase retrieval in the Fresnel region. 2. Partially coherent illumination. Opt Comm. 2006;259:569–580.

23. Guigay J, 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:1617–1619. [PubMed]

24. Wu X, Yan A. Phase Retrieval From One Single Phase Contrast X-Ray Image. Opt Express p Opt Express. 2009;17:11187–11196. [PubMed]

25. Allen L, Oxley M. Phase retrieval from series of images obtained by defocus variation. Opt Comm. 2001;199:65–75.

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

27. Fienup J. Reconstruction of an object from the modulus of its Fourier Transform. Opt Lett. 1978;3:27–29. [PubMed]

28. Fienup J. Phase retrieval algorithms: a comparison. Appl Opt. 1982;21:2758–2769. [PubMed]

29. Dyson N. X-Rays in Atomic and Nuclear Physics. Longman Scientific and Technical; Essex, UK: 1973.

30. Wu X, Dean A, Liu H. Biomedical Photonics Handbook. chap 26. CRC Press; Tampa, Fla: 2003. pp. 26-1–26-34.

31. Hubbell JH, Veigele WI, Briggs EA, et al. Atomic form factors, incohoerent scattering functions, and photon scattering cross sections. Journal of Physical Chemistry Reference Data. 1975;4:471–538.

32. Rudin L. Report #TR:5250:87. Caltech: C, S, Dept; 1987. Images, numerical analysis of singularities and shock filters.

33. Rudin LI, Osher S, Fatemi E. Nonlinear total variation based noise removal algorithms. Physica D. 1992;60:259–268.

34. Wu X, Barnes GT, Tucker DM. Spectral dependence of glandular tissue dose in screen-film mammography. Radiology. 1991;179:143–148. [PubMed]

35. Wu X, Gingold EL, Barnes GT, Tucker DM. Normalized average glandular dose in Molybdenum target-Rhodium filter and Rhodium target-Rhodium filter mammography. Radiology. 1994;193:83–89. [PubMed]

36. Seldin J, Fienup J. Numerical investigation of the uniqueness of phase retrieval. J Opt Soc Am A. 1990;7(3):412–427.

37. Roddier F, Roddier C. Wavefront reconstruction using Iterative Fourier transforms. Appl Opt. 1991;30:1325–1327. [PubMed]

38. Roddier C, Roddier F. Wave-front reconstruction from defocused images and the testing of ground-based optical telescopes. J Opt Soc Am A. 1993;10:2277–2287.

39. Gureyev T, Roberts A, Nugent K. Partially coherent fields, the transport-of-intensity equation, and phase uniqueness. J Opt Soc Am A. 1995;12:1942–1946.

40. Gureyev T, Nugent K. Phase retrieval with the transport-of-intensity equation. II. Orthogonal series solution for nonuniform illumination. J Opt Soc Am A. 1996;13:1670–1682.

41. Tychonoff A, Arsenin V. Solution of Ill-posed Problems. Winston & Sons; Washington: 1977.

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. |