## Abstract

For medical applications of the in-line phase-contrast x-ray imaging, phase retrieval is a crucial step for quantitative imaging such as reconstructing the 3-D distribution of tissue linear attenuation coefficients and refraction indices. The conventional phase retrieval algorithms, such as the transport of intensity equation (TIE) based algorithms, are not robust against the quantum noise that appears in acquired images due to the radiation dose constraints in medical imaging. In this work a new attenuation-partition based iterative phase retrieval algorithm is proposed. This new algorithm takes advantage of the correlation between the attenuation and phase-shift, and is much robust against noise in acquired images. Phase retrieval results from experimental images show that this new iterative algorithm is fast and robust, and it has good potential for medical x-ray imaging applications.

© 2008 Optical Society of America

## 1. Introduction

X-ray imaging is ubiquitously utilized in modern medicine, science and technology. The phasecontrast x-ray imaging is getting more attention recently for its great potentials for medical imaging applications [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]. Phase-contrast x-ray imaging is based on the highly sensitive phase shifts of transmitted x-ray wave field from imaged objects. The amount of phase shift is determined by the refractive index of the tissue. The refractive index n for x-ray is complex and is given by

where *δ*, the refractive index decrement, is responsible for x-ray phase shift, and *β* is the imaginary part of *n* and it is responsible for x-ray attenuation. The *δ* is given by [1, 2, 12]:

where *r _{e}*,

*N*,

_{k}*Z*and

_{k}*f*are the classical electron radius, atomic density, atomic number and real part of the anomalous atomic scattering factor of the element

^{r}_{k}*k*, respectively. Note that tissue’s

*δ*is much larger than

*β*. For example, for the infiltrating ductal carcinoma (IDC) in breast, we found that the cancer tissue’s

*δ*(10

^{-6}–10

^{-8}) is about 1000 times greater than

*β*(10

^{-9}–10

^{-11}) for x-ray of 10 keV to 150 keV [6]. In a word, human tissue’s phase shift-shift cross sections are about 1000 times greater than that of attenuation for x-ray in above energy range [1, 2, 5, 6, 7]. Hence phase-contrast imaging has great potential of improving x-ray imaging sensitivity and tremendously reducing radiation doses to patients for clinical imaging. The in-line phase contrast x-ray imaging is the simplest among all phase contrast x-ray imaging techniques, since no x-ray optics is needed. In this imaging technique, x-ray propagates freely from the imaged object’s exit to the imaging-detector. Phase-shifted x-ray diffracts and generates the dark-bright interference fringes at tissue’s boundary and interfaces [1, 2, 3, 4, 5, 6, 7]. In the inline phase contrast imaging setting, as the x-ray traverses an object, not only its intensity is attenuated, but x-ray undergoes phase-shifts as well. Propagating from the object to the detector a distance downstream, the x-rays with different phase-shifts diffract and interfere to form the dark-bright fringes along the tissue interfaces and edges. It should be noted that x-ray diffraction process from object’s exit to the detector plays the role of x-ray optics to convert x-ray phase-shifts to xray intensity variation for detection. In this way tissue’s boundary and interfaces get enhanced. Furthermore, recovering the phase shift, or the phase retrieval, plays crucial roles in exploring tissue phase contrast. The importance of phase retrieval can be summarized in, and not limited to, the following categories:

- It can disentangle tissue phase-contrast from mixed contrast in a phase-contrast image for phase-contrast enhancement. Tissue’s phase-shift differences are more sensitive than attenuation differences. Enhanced images can be expected based on tissue’s phase-shift compared to that based on attenuation change. This can potentially reduce the x-ray radiation dose, a critical issue for medical imaging fields [5].
- It can provide quantitative maps of tissue projected electron densities. In fact when x-ray traverses an object, the phase shift accrued is given by $\varphi \left(\overrightarrow{\mathbf{r}}\right)=-\frac{2\pi}{\lambda}{r}_{e}\int \delta (\overrightarrow{\mathbf{r}},z)\mathrm{dz},$ , where the integral is over the ray-path. From the formula for
*δ*, Eq. (2), one obtains*ϕ*(**r⃗**)=-*λr*(_{e}∫ρ_{e}**r⃗**,*z*)*dz*=-*λr*(_{e}ρ_{e,p}**r⃗**), where**r⃗**is the position vector at the object plane,*λ*the x-ray wavelength,*ρ*(_{e}**r⃗**,*z*) denotes the object’s electron density,*ρ*(_{e,p}**r⃗**) the integral of tissue electron densities over the ray path, or called the projected electron density as well, and*r*is the classical electron radius._{e} - It is the key step for obtaining 3-D Tomograms of tissue linear attenuation coefficients and refraction indices. This is because that the “direct” phase-contrast tomography without phase-retrieval generates tomograms representing a mixed sum of tissue’s linear attenuation plus a contribution from the 3-D Laplacian of tissue refractive indices, as well as a contribution of artifacts depending on the global distribution of attenuation coefficients and refractive indices. The reconstructed tomographs with this “direct” approach hinders quantitative tissue characterization, and are susceptible to image-artifacts [13]. Therefore, a correct approach should perform the phase retrieval together with filtered backprojection for reconstructing tomograms of tissue linear attenuation coefficients, and tomograms of tissue refraction indices, respectively. Note that tomograms of tissue refraction indices are equivalently tomograms of tissue electron densities [14, 15, 16, 17].

Phase retrieval is to reconstruct the object phase-shift map from multiple phase-contrast images [4, 8, 18, 19]. Phase retrieval is based on x-ray propagation equation derived either from the Fresnel diffraction or the Wigner distribution based phase-space formalism [3, 4, 5, 8, 20]. To be specific, we model the imaged object by a two-dimensional transmission function *T*
_{0}(**r⃗**)=*A*
_{0}(**r⃗**)*e ^{iϕ}*(

**r⃗**), where

*ϕ*(

**r⃗**) denotes the x-ray phase-shift caused by an imaged object, and

*A*

_{0}(

**r⃗**) denotes the x-ray transmission, or the attenuation-map of the object. For the imaging geometry, we assume that the source-object distance is

*R*

_{1}, and the object-detector distance is

*R*

_{2}, hence the magnification factor for the object image is

*M*=(

*R*

_{1}+

*R*

_{2})/

*R*

_{1}. We stress again no x-ray-optics such as gratings or zone-plates are used in imaging. A general formula for x-ray intensity at the detector was derived as [5]:

$$\phantom{\rule{.9em}{0ex}}\phantom{\rule{.5em}{0ex}}+\left[2\mathrm{sin}\left(\frac{\pi \lambda {R}_{2}}{M}{\overrightarrow{\mathbf{u}}}^{2}\right)-\left(\frac{2\pi \lambda {R}_{2}}{M}{\overrightarrow{\mathbf{u}}}^{2}\right)\xb7\mathrm{cos}\left(\frac{\pi \lambda {R}_{2}}{M}{\overrightarrow{\mathbf{u}}}^{2}\right)\right]\xb7\hat{\U0001d4d5}\left({A}_{0}^{2}\varphi \right)-$$

$$\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}\phantom{\rule{.2em}{0ex}}-\mathrm{cos}\left(\frac{\pi \lambda {R}_{2}}{M}{\overrightarrow{\mathbf{u}}}^{2}\right)\xb7\frac{\lambda {R}_{2}}{2\pi M}\xb7\hat{\U0001d4d5}\left(\nabla \xb7\left({A}_{0}^{2}\nabla \varphi \right)\right)\},$$

where,

*I*
_{in} is the incident x-ray intensity at *R*
_{1}, *λ* the wavelength of the monochromatic point x-ray source and $\widehat{\U0001d4d5}$
(·) is the Fourier transform. Here *πλR*
_{2}
**u⃗**
^{2}/*M* is the FT-Space Fresnel propagator’s phase. When the FT-Space propagator’s phase is small

Eq. (5) imposes the moderate variation on the object and imaging geometry. In other words, Eq. (5) imposes conditions such that the imaged object’s phase variation should be moderate or the distance from the object to detector cannot be too large. In clinical imaging, the spatial resolution is limited (<20–30 mm/lp) and object-detector distance <3 m, so the condition Eq. (5) is satisfied. But for synchrotron-based phase contrast imaging experiments the Eq. (5) can be violated [20].

Equation (3) can be simplified to the common known Transport of Intensity Equation (TIE), [21, 18, 5]

So, TIE (6) is a special case of Eq. (3). It is valid only for moderate resolution images. For high resolution images, i.e. (5) not satisfied, one has to use (3) for phase retrieval. From Eq.(3) or (6), multiple but at least two images (such as *A*
^{2}
_{0} and *I*) should be known. For clinical applications the phase retrieval based on just two acquired images are much preferable than that based on multiple images because of the speed and radiation dose consideration. There are two ways to acquire two images of the imaged object. One way is to make two exposures with different object-detector distances. Another way is to make just one exposure for two images with two detectors in downstream, and we called this approach the dual detector approach [8].

Note (6) is a second order elliptic partial differential equation about phase *ϕ*. The solution *ϕ*, represented in the form of pseudo-differential operator, can be approximated by (see [19] for details about the implementation of the algorithm)

where ∇ is the gradient operator. Since ∇^{-2} is singular at zero-frequencies, the algorithm based on (5) is unstable near zero-frequencies with noisy data, since the low frequency noise components will get amplified in this case. In medical imaging applications the noise embedded in images is intrinsic because of the radiation dose constraints. Therefore, the algorithms based on Eq. (7) is not suitable for clinical imaging applications [22]. To overcome this difficulty, some kind of regularization can be used. The most widely used regularization is Tikhonov regularization:

where *α* is called the Tikhonov regularization parameter [23, 24, 25, 22]. In essence Tikhonov regularization seeks the minimum-norm, least squares solution of Eq. (7). However the performance of Tikhonov regularization is highly dependent on the level of noise. Note that in Eq. (6), phase *ϕ* appears only in the term contains *λR*
_{2}/2*πM*. Since the x-ray’s short wavelength *λ*, and clinical restriction (small *R*
_{2}), the contribution of tissue’s phase *ϕ* to intensity is small. A small amount of noise compared to intensity will render Tikhonov regularization useless for phase retrieval [22]. Another disadvantage of Tikhonov regularization is that *α* is not *a priori* parameter which makes it difficult for practical applications.

As we mentioned earlier that the x-ray phase-shift by tissues arises from the x-ray coherent scattering by tissues. In contrast x-ray attenuation arises from three x-ray-tissue interactions: the photoelectric absorption, coherent scattering and incoherent scattering in the x-ray energy range employed for medical imaging. For high energy x-rays of 60 keV to 500 keV x-ray attenuation by soft tissues is determined by the incoherent x-ray scattering. Therefore, we observed that both the x-ray phase-shift and attenuation by soft tissues are all determined by the projected electron density for these high energy x-rays. We called this complementary relationship between phase-shift and attenuation for soft tissues as the phase-attenuation duality.

Looking for more effective and robust phase retrieval algorithms, we found a phase retrieval approach based on just a single acquired phase-contrast image. This approach, which is based on the phase-attenuation duality, is applicable for soft tissue (or light element objects) imaging with x-rays of 60 keV or higher [12]. This phase retrieval approach is very much robust against image noise [22]. In this work, we extend the notion of phase-attenuation duality for light element objects to the correlation between phase shift and attenuation for general objects and general x-ray energy range. Based on this observation of the correlation between phase shift and attenuation, we develop an iterative phase retrieval algorithm, called Attenuation-Partition Based Algorithm (APBA) for general objects. Before we are going to detailed analysis in following sections, we briefly describe how this new algorithm works. Figure 1 represents the flow chart of the algorithm. In the chart, *I*
_{0} and *A*
_{0} represent the acquired phase-contrast image intensity at the source to imaging-detector distance *SID*=*R*
_{1}+*R*
_{2} and the attenuation data from the acquired image at *SID*=*R*
_{1}. The new algorithm consists of the following four simple steps: (1) compute *A*
_{n+1} and *ϕ*
_{n+1} by assuming *I _{n}* satisfy the phase-attenuation duality Eq. (4) of [12]; (2) compute the residue

*δ*

*A*

_{n+1}=

*A*

_{n+1}-

*A*

_{0}; (3) propagate $\delta {A}_{n+1}\xb7{e}^{i{\varphi}_{n+1}}$ using Fresnel propagation to get

*δ*

*; (4) replace*

*I*_{n+1}*I*with ${I}_{n+1}={\left(\sqrt{{I}_{0}}+\sqrt{\delta {I}_{n+1}}\right)}^{2}$ . The phase retrieval is complete through about 10 iteration steps as is shown below.

_{n}## 2. The Attenuation-Partition Based Algorithm (APBA)

In the diagnostic x-ray imaging, the tissue’s attenuation change arises from three x-ray–tissue interactions: the photoelectric absorption, the coherent scattering, and the incoherent scattering [26, 27, 5, 12]. If we denote the x-ray wave amplitude attenuated by incoherent scattering by *A*
_{incoh}, and the amplitude attenuated by photoelectric absorption and coherent scattering by *A*
_{pe,coh}, the total attenuated x-ray amplitude *A*
_{0} can be written as

The reason for factoring out the incoherent scattering contribution is that both *A*
_{incoh} and the phase-shift *ϕ* are all determined by tissue’s electron density. In fact, 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. We noted that in these cases attenuation and phase are both determined by object’s electron densities, we call this notion as the phase-attenuation duality [12]. With the help of the duality notion, we found a phase retrieval approach based on just a single acquired phase-contrast image [12]. Surprisingly this phase retrieval approach is very much robust against image noise [22]. For imaging tasks with general objects containing high z-elements such as iodinated-contrast media, or imaging tasks with x-rays of lower energies such as in mammography, the duality condition will not hold any more, but the correlation between phase shift and attenuation for general objects and general x-ray energy range can be used for phase retrieval through the attenuation-partition as is shown in follows.

First let’s denote the amplitude attenuated by Compton scattering only as *A*
_{KN}. As we pointed out in [12], *A*
_{KN} has a relation with tissue’s phase-shift, *ϕ*, in the form:

where *µ _{KN}* is the tissue linear attenuation coefficient from the Compton scattering,

*λ*the x-ray wavelength,

*r*the classic electron radius and

_{e}*σ*

_{KN}is the total cross section for x-ray photon Compton scattering from a single free electron, which can be determined by the time-honored Klein-Nishina formula:

Here *η*=*E*
_{photon}/*m _{e}c*

^{2}, with

*E*

_{photon}denotes the photon energy of the primary x-ray beam and

*m*

_{e}c^{2}the resting electron energy. When the phase-attenuation duality holds, combining (10) and (6), we get the duality reversion formula [12]:

where *I* is the intensity of phase-contrast image acquired at *SID*=*R*
_{1}+*R*
_{2},

and $\widehat{\U0001d4d5}$ (·), $\stackrel{\u02c7}{\U0001d4d5}$ (·) represent the Fourier and Inverse Fourier transforms. And phase can be retrieved from equations (10) and (12). For imaging tasks with general objects and using x-rays of lower energies, Eq. (10) will not hold any more. In order still to take advantage of the correlation between phase shift and attenuation, we take a strategy of attenuation-partitioning and iterative improvement of approximation. Specifically we rewrite Eq. (9) in the form

where

Now we construct an iterative algorithm as follows. (1) First, we regard the acquired phasecontrast image as if it would receive all contribution from the Compton scattering. In other words, the amplitude √*I* of the phase-contrast image is contributed all from *Ã*(**r⃗**)=*A*
_{KN}(**r⃗**). Using the duality inversion formula (12) compute *Ã* and the approximated phase-map
$\varphi =\frac{\lambda {r}_{e}}{{\sigma}_{\mathrm{KN}}}\mathrm{ln}\left({\stackrel{~}{A}}^{2}\right)$
; (2) second, we compute the error in estimation of the object’s attenuation map from this approximation. This error is apparently *δA*=*Ã*-*A*
_{0}, where *A*
_{0} is the object’s attenuation obtained from the image acquired at *SID*=*R*
_{1}; (3) third, we want to improve our previous estimate on *Ã*(**r⃗**)=*A*
_{KN}(**r⃗**) by incorporating the error *δA*. To do so we consider a hypothetical object *δ*
*T*
_{0}(**r⃗**)=*δA*(**r⃗**)·*e*
^{iϕ(r⃗)}, and compute the transmitted x-ray wave field arriving at the detector at *SID*=*R*
_{1}+*R*
_{2}. The computation is according to the paraxial Fresnel diffraction theory, for a monochromatic x-ray point source of wavelength *λ* at the origin, the diffracted spherical x-ray wave field arriving at **r⃗** on the detector plane, can be modelled by the *Fresnel-Kirchhoff* integral [5]

where
$i=\sqrt{-1}$
,**r⃗** is the position in the imaged object, *R*
_{1} the distance from point x-ray source to object, *R*
_{2} the distance from object to detector plane, *I*
_{in} the incident x-ray intensity at *R*
_{1}. For convenience we denote

$${E}_{1}\left(\overrightarrow{\mathbf{r}}\right)=\frac{M}{\sqrt{{I}_{\mathrm{in}}}}E(M\overrightarrow{\mathbf{r}};{R}_{1},{R}_{2})\xb7{e}^{-\frac{i}{2\omega}h\left(\overrightarrow{\mathbf{r}}\right)},$$

where *M*=(*R*
_{1}+*R*
_{2})/*R*
_{1} is the geometric magnification factor. Note in above *E*
_{1}(**r⃗**) is the normalized x-ray wave field at detector entrance, and its intensity forms the phase-contrast image. From Eq. (15) a simple computation leads to

with *G*(**r⃗**)=exp[*i*
**r⃗**
^{2}/2*ω*],∗ represents the convolution operator, and
$\mid {E}_{1}\mid =\frac{M}{\sqrt{{I}_{\mathrm{in}}}}\sqrt{I}$
. In the rest of the paper,
$\frac{{M}^{2}}{{I}_{\mathrm{in}}}I$
will be called the normalized intensity and still be denoted by *I*. Rewrite
${E}_{1}\left(\overrightarrow{\mathbf{r}}\right)=\sqrt{I\left(\overrightarrow{\mathbf{r}}\right)}{e}^{i\psi \left(\overrightarrow{\mathbf{r}}\right)}$
, we found by applying stationary phase formula [28], that

Eq. (17) provides a general formula for the x-ray wave eikonal transport. Now for given object phase *ϕ*, denoting two virtual objects, *T̃*
_{0}=*Ã*
*e ^{iϕ}* and

*δ*

*T*

_{0}=

*δA*located at the object plane. After the paraxial Fresnel free propagation to position

^{eiϕ}*SID*=

*R*

_{1}+

*R*

_{2}, i.e. Eq. (16), the corresponding propagated wave fields at the detector-entrance are denoted by

*Ẽ*

_{1}=√

*Ĩe*, and $\delta {E}_{1}=\sqrt{\delta I}{e}^{i\delta \psi}$ respectively. Then by (16),

^{i$\tilde{\psi}$}*E*

_{1}=

*Ẽ*

_{1}-

*δ*

*E*

_{1}. This gives

From (17)

$\stackrel{~}{\psi}\left(\overrightarrow{\mathbf{r}}\right)-\delta \psi \left(\overrightarrow{\mathbf{r}}\right)=\frac{\omega}{2}\left(\frac{{\nabla}^{2}\stackrel{~}{A}\left(\overrightarrow{\mathbf{r}}\right)}{\stackrel{~}{A}\left(\overrightarrow{\mathbf{r}}\right)}-\frac{{\nabla}^{2}\delta A\left(\overrightarrow{\mathbf{r}}\right)}{\delta A\left(\overrightarrow{\mathbf{r}}\right)}\right)$ .

Thus

$\mathrm{cos}\left(\stackrel{~}{\psi}\left(\overrightarrow{\mathbf{r}}\right)-\delta \psi \left(\overrightarrow{\mathbf{r}}\right)\right)=\mathrm{cos}\left(\frac{\omega}{2}\left(\frac{{\nabla}^{2}\stackrel{~}{A}\left(\overrightarrow{\mathbf{r}}\right)}{\stackrel{~}{A}\left(\overrightarrow{\mathbf{r}}\right)}-\frac{{\nabla}^{2}\delta A\left(\overrightarrow{\mathbf{r}}\right)}{\delta A\left(\overrightarrow{\mathbf{r}}\right)}\right)\right)\approx 1$

when we assume

${\left[\frac{\omega}{2}\left(\frac{{\nabla}^{2}\stackrel{~}{A}\left(\overrightarrow{\mathbf{r}}\right)}{\stackrel{~}{A}\left(\overrightarrow{\mathbf{r}}\right)}-\frac{{\nabla}^{2}\delta A\left(\overrightarrow{\mathbf{r}}\right)}{\delta A\left(\overrightarrow{\mathbf{r}}\right)}\right)\right]}^{2}\ll 1$ .

This gives

by assuming *Ĩ*(**r⃗**)>*δI*(**r⃗**). In other words, the square root of the (normalized) phase-contrast image intensity data *I*, generated by the object, is the square root of the intensity *Ĩ*, generated by the virtual object *T̃*
_{0}=*Ãe ^{iϕ}*, minus the square root of the intensity

*δI*, generated by the virtual object

*δT*

_{0}=

*δAe*. This elaborated analysis leads to the 4th step of the algorithm, that is, (4). update √

_{iϕ}*Ĩ*with $\sqrt{I}+\sqrt{\delta I}$ . A new iteration step can then be started. Summarizing the above analysis, we arrive at the following new phase-retrieval algorithm:

**Algorithm 1**. *The Attenuation-Partition Based iterative phase retrieval Algorithm (APBA). Assume the attenuation A _{0} and the normalized intensity I*=

*M*

^{2}/

*I*

_{in}·|

*E*(

**r⃗**;

*R*

_{1},

*R*

_{2})|

^{2}

*are known inputs. The initial value of δI, generally 0, is also set from the input.*

*1. Compute Ĩ from (19).*

*2. Compute*
*Ã*=A_{KN}
*from Eq. (12) with*
$\frac{{M}^{2}}{{I}_{\mathrm{in}}}$
*I in (12) substituted with Ĩ. Then compute ϕ from (10).*

*3. Compute*
*δA*=*Ã*-*A*
_{0}.

*4. Compute the improvement of the iteration:*

*where ‖·‖ denotes some kind of norm, that can effectively reflect the accuracy of the retrieved data as well as the visual perception of the data as an image, and*
*δA*
_{old}
*is the δA computed in the previous iteration step. If δ is less than a predefined threshhold value β, or the iteration step exceeds a predefined maximum number of iteration steps, then ϕ is the retrieved phase and the iteration stops; otherwise, go to the next step.*

*5. Compute δI=|E _{1}|^{2} from (16), with T*

_{0}=

*δA·e*.

^{iϕ}and go back to step (1)Remarks:

- Due to space limitation a detailed mathematical proof about this algorithm will be present in another paper. It can be shown that this algorithm converges whenever the parameter
*k*in (12) is greater than zero. - The extent of correlation between
*A*_{0}and*ϕ*is another important factor affecting the speed of convergence. In other words, the relative portion of incoherent scattering contribution in the total attenuation can affect the speed of convergence. In real situations, if*k*is chosen correctly, an accurate result can be reached within 10 iteration steps.

## 3. Phase retrieval with experimental images

In this section, we compare the performance of APBA with the TIE algorithm [19] for lab experimental images. Two groups of data are tested. In the first group, the two images are that of an acrylic plate of 1.42 mm. Acrylic is a clear plastic of C_{5}H_{8}O_{2} that resembles glass, and is with commercial names such as Lucite and Plexi glass. Its mass density ranges from 1.16–1.20 g/cm3. The laser-cut sharp edge of the acrylic plate was imaged with a micro-focus x-ray tube (UltraBright Microfocus X-ray source, Oxford Instruments, CA) operating at 40 kVp. The focal spot used was of 14 *µm*. For the imaging detector a Computed Radiography(CR)-system was used, which includes a CR-plate (the photostimulable phosphor) and a reader of a 43.75 *µm* sampling pitch (Konica REGIUS 190, Konica Minolta, Japan). Figure 2(a) and Fig. 2(b) are the central parts of the acquired images. Figure 2(a) is a conventional radiograph of the plate acquired with a source-detector distance *SID*=66 cm, and an object-detector distance *R*
_{2}≈0. Figure 2(b) is a phase-contrast image acquired with *SID*=132 cm, and an object-detector distance *R*
_{2}=66 cm. In Fig. 2(b) the phase-contrast manifests itself as the dark-bright fringe along the acrylic plate-edge. The dark-bright fringe pair is the signature enhancement-pattern of phase-jumps. It forms as a result of the interference of the diffracted x-ray with different phase-shifts along the plate edge. As a comparison, the corresponding image acquired with *R*
_{2}≈0, that is, the conventional radiograph of the plate (Fig. 2(a)), showed no dark-bright fringe. This is because x-rays exiting from the plate need to diffract a certain distance to generate interference between waves of different phase-shifts [6, 10]. We cut a piece with 1500×1500 pixels from the original image data (*SID*=132 cm, *R*
_{2}=66 cm) and treat it as the intensity data *I*. Based on the magnification factor *M*, we cut a corresponding part from the image with (*SID*=66 cm, *R*
_{2}≈0) and treat it as the attenuation data *A*
^{2}, and interpolate it to the same size as *I*. Care is taken in the process of aligning the two images. Figure 2(c) is the recovered phase map using APBA after 10 iteration steps. The corresponding phase maps recovered using TIE algorithm [19], with different Tikhonov regularization parameter *α*, are shown in Fig. 3. We can see that the result using APBA is much better, in visual perception, than that using TIE algorithm [19]. The true phase shifts of the acrylic plate is about 101 radian, which is estimated from the estimated average x-ray wavelength, the plate’s composition and mass density for acrylic, and the plate thickness (1.5 mm). We found the retrieved phase shifts for 10 iteration steps is 57.3, about half the estimated “intrinsic” phase-shift. The error is attributed to the image alignment error, and the error in mass-density estimation of the plate. On the other hand, the estimated phase-shifts of the plate from the TIE-based retrieval algorithm (Fig. 3) are 22941, 36.1 and 2.6 when *α*=0, 10Δ^{2}, and 50Δ^{2} respectively, where Δ is the sampling step-size in FT-space. Hence it is clear that the phase-maps by TIE based retrieval are way off the true phase-shift values.

In the second group, two images of a breast cancer patient’s lumpectomy specimen were used for testing the phase retrieval algorithms. As the details of image acquisition were described in [29], here we just give a brief description. The source used was a micro-focus x-ray tube (Model L8121-01, Hamamastu Photonics, Japan) operating at 40 kVp. The focus spot was set to 7*µm*. The detector is an amorphous-selenium based flat panel detector (DirectRay, Hologic, DE) with a pixel pitch of 139*µm*. Figure 4(a) is the conventional radiograph of the lumpectomy specimen acquired with *SID*=66 cm, *R*
_{2}≈0. Figure 4(b) is the phase-contrast image of the specimen acquired with *SID*=184.9 cm and *R*
_{2}=118.5 cm. We used these two images as the input for testing the phase retrieval algorithms. Compared to Fig. 4(a), Fig. 4(b) shows the enhanced dark-bright fringes at tissue boundaries. The recovered phase map using APBA (after 10 iteration steps) is shown in Fig. 4(c). However the phase recovery using TIE based algorithm, as shown in Fig. 5, was completely failed. This case shows again our new algorithm is much more robust than TIE algorithm.

## 4. Conclusions

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. The conventional phase retrieval algorithms, such as the transport of intensity equation (TIE) based algorithms, are not robust against the quantum noise that appears in acquired images. In this work a new attenuation-partition based iterative phase retrieval algorithm is proposed. This new iterative algorithm is derived from our previous notion of the phase-attenuation duality, and it takes advantage of the correlation between the attenuation and phase-shift. In this work, we analyzed the speed of convergence of this algorithm and conclude that the converge speed is proportional to the order of the correlation between the objects attenuation and phase. The phase retrieval with the experimental images showed that the new algorithm has much better performance than the TIE-based algorithm ([21]). We expect that this new iterative phase retrieval will find many applications in medical x-ray imaging.

## Acknowledgments

The authors would like to thank Mr. Da Zhang for his help with the experimental images.

## References and links

**1. **A. Snigirev, I. Snigireva, V. Kohn, S. Kuznetsov, and I. Shelokov, “On the possibilities of x-ray phase contrast microimaging by coherent high-energy synchrotron radiation,” Rev. Sci. Instrum. **66**, 5486–5492 (1995). [CrossRef]

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

**3. **A. Pogany, D. Gao, and S. Wilkins, “Contrast and resolution in imaging with a microfocus x-ray source,” Rev. Sci. Instrum. **68**, 2774–2882 (1997). [CrossRef]

**4. **D. Paganin and K. Nugent, “Noninterferometric Phase Imaging with Partially Coherent Light,” Phy. Rev. Lett. **80**, 2586–2589 (1998). [CrossRef]

**5. **X. Wu and H. Liu, “A general theoretical formalism for X-ray phase contrast imaging,” J. X-ray Sci. Technol. **11**, 33–42 (2003).

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

**7. **E. Donnelly, R. Price, and D. Pickens, “Dual focal-spot imaging for phase extraction in phase-contrast radiography,” Med. Phys. **30**, 2292–2296 (2003). [CrossRef] [PubMed]

**8. **X. Wu and H. Liu, “A dual detector approach for X-ray attenuation and phase imaging,” J. X-ray Sci. Technol. **12**, 35–42 (2004).

**9. **X. Wu and H. Liu, “An experimental method of determining relative phase-contrast factor for x-ray imaging systems,” Med. Phys. **31**, 997–1002 (2004). [CrossRef]

**10. **X. Wu and H. Liu, “A new theory of phase-contrast x-ray imaging based on Wigner distributions,” Med. Phys. **31**, 2378–2384 (2004). [CrossRef] [PubMed]

**11. **X. Wu and H. Liu, “A reconstruction formula for soft tissue X-ray phase tomography,” J. X-ray Sci. Technol. **12**, 273–279 (2004).

**12. **X. Wu, H. Liu, and A. Yan, “X-ray phase-attenuation duality and phase retrieval,” Opt. Lett. **30**(4), 379–381 (2005). [CrossRef] [PubMed]

**13. **X. Wu, H. Liu, and A. Yan, “Phase-Contrast X-Ray Tomography: Contrast Mechanism and Roles of Phase Retrieval,” Eur. J. Radiol. (to be published)(2008). [CrossRef] [PubMed]

**14. **A. Bronnikov, “Theory of quantitative phase-contrast computed tomography,” J. Opt. Soc. Am. A **19**, 472–480 (2002). [CrossRef]

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

**16. **T. Gureyev, D. Paganin, G. Myers, Y. Nesterets, and S. Wilkins, “Phase-and-amplitude computer tomography,” Appl. Phys. Lett. **89**, 034,102 (2006). [CrossRef]

**17. **R. M. P. Cloetens, M. Schlenker, and S. Lerbs-Mache, “Quantitative phase tomography of Arabidopsis seeds reveals intercellular void network,” PNAS **103**, 14,626–14,630 (2006). [CrossRef]

**18. **K. Nugent, T. Gureyev, D. Cookson, D. Paganin, and Z. Barnea, “Quantitative Phase Imaging Using Hard X Rays,” Phy. Rev. Lett. **77**, 2961–2964 (1996). [CrossRef]

**19. **L. Allen and M. Oxley, “Phase retrieval from series of images obtainedbe defocus variation,” Opt. Commun. **199**, 65–75 (2001). [CrossRef]

**20. **X. Wu and H. Liu, “Phase-space evolution of x-ray coherence in phase-sensitive imaging,” Appl. Opt. **47**, E44–E52 (2008). [CrossRef] [PubMed]

**21. **M. Teague, “Deterministic phase retrieval: a Green’s function solution,” J. Opt. Soc. Am. **73**, 1434–1441 (1983). [CrossRef]

**22. **X. Wu and H. Liu, “Robustness of a phase-retrieval approach based on phase-attenuation duality,” J. X-ray Sci. and Tech. **15**, 85–95 (2007).

**23. **A. Tikhonov and V. Arsenin, *Solution of Ill-posed Problems* (Winston & Sons, Washington, 1977).

**24. **I. Schelokov, T. Weitkamp, and A. Snigirev, “Reconstruction of an object phase transmission function from inline X-ray holograms,” Opt. Commun. **213**, 247–258 (2002). [CrossRef]

**25. **T. Gureyev, A. Pogany, D. Paganin, and S. Wilkins, “Linear algorithms for phase retrieval in the Fresnel region,” Opt. Commun. **231**, 53–70 (2004). [CrossRef]

**26. **N. Dyson, *X-Rays in Atomic and Nuclear Physics* (Longman Scientific and Technical, Essex, UK, 1973).

**27. **X. Wu, A. Dean, and H. Liu, *Biomedical Photonics Handbook*, chap. 26, pp. 26-1–26-34 (CRC Press, Tampa, Fla., 2003).

**28. **L. Hörmander, *The Analysis of Linear Partial Differential Operators I* (Springer, New York, 1983). [CrossRef]

**29. **D. Zhang, M. Donvan, L. Fajardo, A. Archet, X. Wu, and H. Liu, “Preliminary Feasibility study of an in-line phase contrast x-ray imaging prototype,” IEEE Trans Biomedical Engineering **55**, (in press) (2008).