Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Opt Express. Author manuscript; available in PMC 2010 August 23.
Published in final edited form as:
PMCID: PMC2925321

Performance analysis of the attenuation-partition based iterative phase retrieval algorithm for in-line phase-contrast imaging


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.

1. Introduction

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 [var phi] by tissues is determined by


where re 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 [var phi] (r) be the x-ray phase-shift caused by the imaged object, and A0(r) the x-ray transmission, or the attenuation-map of the object. Then the Fourier transform of the x-ray intensity, F (I)(u) = R2 I(x)exp [2 πix · u]dx, at position away from the object with distance R2, of the monochromatic point x-ray source starting at a place away from the object R1 distance, can be modeled by the following [9]

equation image

where Iin is the incident x-ray intensity at R1, λ the wavelength of the monochromatic point x-ray source and M = (R1 + R2)/R1 is the geometric magnification. When the FT-Space Fresnel propagator πλ R2u2/M [double less-than sign] 1, Eq. (2) can be simplified to the Transport of Intensity Equation (TIE) [21, 4, 9]


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 πλ R2u2/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.

2. The attenuation-partition based algorithm (APBA) and an image accuracy measure

2.1. The Attenuation-Partition Based Algorithm

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 A0 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 AKN, which is dominated by the x-ray Compton scattering, deserves a special attention. This is because both of AKN and the x-ray phase-shift [var phi] are determined by the tissue’s projected electron density:


where λ is the x-ray wavelength, re 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:


with η = Ephoton/mec2. Here we denote the photon energy of the primary x-ray beam by Ephoton and mec2 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. A0AKN[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 R1 from the source. We denote R2 as the distance from object to detector plane, M = (R1 + R2)/R1 the geometric magnification factor, Iin the entrance x-ray intensity at R1, and ID(rD) the x-ray intensity at the detector plane. For convenience, we denote I = M2 · ID(rD)/Iin as the normalized intensity of phase-contrast image. When the phase-attenuation duality holds, the phase map [var phi] (r) can be robustly retrieved from just a single projection image[14]:

equation image



and An external file that holds a picture, illustration, etc.
Object name is nihms226357ig1.jpg, 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 A0 as


where we denote the attenuation caused by photoelectric absorption and coherent scattering by Ape,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 A0 = AKN · Ape,coh (Eq. (8)) we actually factor the small residual binding effect of atomic electrons into Ape,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 A02 acquired with R2 = 0, and the other is the acquired phase contrast image I = M2 · ID(rD)/Iin with R2 > 0. Employing the acquired phase contrast image I and the duality transform Eq. (6), we will first obtain an estimate for the attenuation-component AKN and phase map [var phi]. We then rewrite Eq. (8) as


and find the correction term δ A using the estimate of AKN. We then employing the Fresnel Diffraction transform (as defined in Eq. (11)) to transport the wave field δ Aei[var phi] from R1 to R2. We can find δI = |An external file that holds a picture, illustration, etc.
Object name is nihms226357ig2.jpg (δAei[var phi])|2, which is the difference between phase contrast image I and the “duality-only” counterpart IKN = |An external file that holds a picture, illustration, etc.
Object name is nihms226357ig2.jpg (AKNei[var phi])|2. With the corrected “duality-only” image intensity IKN=(I+δI)2 we can start a new round of iterations by repeating above procedure. For a rigorous analysis of the iterative algorithm and its convergence interesting readers are referred to [1]. Note that the equation I=IKNδI is generally valid, since it is actually a result of x-ray Fresnel diffraction and extremely smallness of hard x-ray wavelength compared to finest resolution achievable in the phase contrast imaging. While interesting readers can find the mathematical proof of this equation in [1], an intuitive explanation of this formula comes from the x-ray propagation. In such a wave propagation process, or the so-called semiclassical wave propagation, the phase of a wave field evolves essentially according to the free-space Hamilton-Jacobi equation from its initial phase value. So if we denote the Fresnel free propagation as a Fresnel transform An external file that holds a picture, illustration, etc.
Object name is nihms226357ig2.jpg acting on the initial wavefield, therefore the two resulted wavefields An external file that holds a picture, illustration, etc.
Object name is nihms226357ig2.jpg (AKN exp[i[var phi]]) and An external file that holds a picture, illustration, etc.
Object name is nihms226357ig2.jpg (δ A exp[i[var phi]]) would have the same resultant phases, so the resultant amplitude from the two-wave superposition is given as I=IKNδI.

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

Fig. 1
Flow chart of APBA


In an imaging experiment, assume we have performed two imaging measurements. One image (radiograph) is the attenuation image A02 acquired at SID= R1, another is a normalized phase-contrast image I = M2 · ID(rD)/Iin acquired at SID= R1 + R2. With A02 and I as well as the initial δ I, usually 0, as known inputs, we first assume IKN=(I+δI)2. Then

  1. Compute AKN=D(IKN) and [var phi] from the duality transform Eq. (6).
  2. Compute δ A from
    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 = Ape,coh in the ideal case and Ape,coh is bounded between 1 and A0. 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(A0), the minimum value of A0, to P in the iterative computations, we can make the algorithm more stable. Moreover, if we know a better lbd for Ape,coh, other than the minimum of A0, the converging speed of the algorithm can be greatly improved.
  3. Compute δI by Fresnel propagate δAei[var phi] from position R1 to R2: δI = |An external file that holds a picture, illustration, etc.
Object name is nihms226357ig2.jpg (δAei[var phi]|2 with
  4. Compute IKN=(I+δI)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 ||δIk+1δIk|| is less than a predefined threshold value β, or the iteration step exceeds a predefined maximum number of iteration steps, then [var phi] is the retrieved phase and the iteration stops; otherwise further iteration is needed.

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

2.2. An Image Accuracy Measure

A continuous signal is generally represented as a function of vector variables: f(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 (r), 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:


where Ω is the finite domain of image functions f and g, V (Ω) represents the area of domain Ω, and μ = Ω(gf)dr/V (Ω) is the statistical mean value of gf. 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 L1 norm of the gradient of the image function, and not the L2 norm. For two image functions f and g, we define the TV norm of gf as the closeness of the two images:


where [nabla] is the gradient operator. Since std is the form of L2 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 [var phi]1[var phi]2 and [var phi]1[var phi]3 are std([var phi]1, [var phi]2) = 0.206, std([var phi]1, [var phi]3) = 0.472 respectively. But the image [var phi]3 is much more closer to [var phi]1 than [var phi]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 [var phi]1[var phi]2 and [var phi]1[var phi]3 are TV([var phi]1, [var phi]2) = 0.0247 and TV([var phi]1, [var phi]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.

Fig. 2
Example “Lena” images used in measuring the closeness between images. (a) [var phi]1, (b) [var phi]2, (c) [var phi]3

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 [var phi] is unique up to a constant with given information about the attenuation-map A0 and the phase-contrast intensity-map I.

3. Simulation Tests

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, R1, and object to detector, R2, are set to 26 inches (0.66 m) respectively. In this way the magnification factor M = (R1 + R2)/R1 is equal to 2. For convenience, the incident x-ray intensity Iin at R1 is set to M2 (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).

3.1. A breast specimen model

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 Ape,coh2 and AKN2 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.

Fig. 5
Image representation of the inputs generated from the simulation model and Fresnel propagation. (a) the phase map [var phi]; (b) the attenuation map A02; and (c) the normalized Fresnel propagated phase contrast image I with object to detector distance ...

The attenuation image A02 of the specimen model simulated with 35.5 keV x-ray, and its two corresponding partition images AKN2 and Ape,coh2, are shown in Fig. 3. For a comparison, the profiles of AKN2 and Ape,coh2 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 Ape,coh2 to the total attenuation A02 gets smaller when x-ray energy is getting higher. Especially, when x-ray energy equals 59.5 keV, the contribution of Ape,coh2 for the soft tissue can almost be neglect, as is expected.

Fig. 3
Image manifest of (a) Ape,coh2, (b) AKN2 and (c) A02 when x-ray energy equals 35.5 keV.
Fig. 4
Profiles of Ape,coh2, the solid lines, and AKN2, the dotted lines, when x-ray energy equals (a) 18.5 keV, (b) 35.5 keV and (c) 59.5 keV.

3.2. Comparison with the GS Algorithm

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 R2 ≈ 1 m, and converges faster for larger R2, 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 R1 = R2 = 26 inches (0.66 m) respectively. For convenience, the incident x-ray intensity Iin at R1 is set to M2, where M = (R1 + R2)/R1 is the magnification factor. The phase map [var phi] and the attenuation A02 are generated from our phantom model for 35.5 keV x-ray. Fig. 5 shows the simulated images of the phase map [var phi], the attenuation image A02 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.04E – 3, almost 33 times greater than that, 3.17E – 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.

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

3.3. Comparison with the Transport of Intensity (TIE) algorithm

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 [var phi] is given by


for the given normalized phase-contrast image I and attenuation image A02. Here [nabla] is the gradient operator, and [nabla]−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.

Fig. 7
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 [nabla]−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 [nabla]−2 by |u|2/(|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 [var phi]κ corresponding this κ is as close to the true phase [var phi] as possible. Roughly speaking, the regularization parameter κ is inversely proportional to the images signal-noise ratio. Two problems, however, arise to this regularization. First, the true phase [var phi] is not known in real situations. So the regularization parameter κ is not a priori, which makes it difficult in practical applications. Second, this Tikhonov regularization is based on finding a stable solution to Ax = y, in Hilbert spaces X, Y, by solving the minimization problem


It is L2 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 A02 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 = σb/Nb corresponding to the direct exposure area (where A2 = 1) outside the object in background. The Poisson statistics dictates that Nb=1/δb2 and in this way the photon count N(i, j) can be determined accordingly at each pixel. Once N(i, j) is determined, the statistical errors at each pixel can be assigned using a computer simulated random Poisson distribution generator with mean corresponding to the photon counts N(i, j). In the simulations below, the background noise level δb is set to 0.0003. The images with the noise added are shown in Fig. 8(b) and (c). The quality of an image depends not only on the noise level but also on the extent of image contrast change. With the assumption P(i, j) = c · N(i, j) above, one can easily see that δb=δPb/Pb¯ is the statistical coefficient of variation in absorbed photon numbers in background. δoi=δPoi/Poi¯, the coefficient of variation of the object image, on the other hand, is the structural coefficient of variation of the sampled image, which reflects the image’s normalized extent of image contrast change. We will use the ratio δoi/δb to reflect the image quality. The larger the ratio, the higher the image quality. In our simulation models, when δb = 0.0003, the corresponding ratio δoi/δb for attenuation A02 and phase-contrast image I, Fig. 8(b) and (c), are 1.67 and 6.35 respectively. Because of the phase-contrast effect, the image quality of the phase-contrast image I is higher than the attenuation image A02 although they have the same background noise level.

Fig. 8
Comparison of the TIE algorithm and APBA with noise added. (a) True phase map [var phi]; (b) attenuation map A02; (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 A02 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.

Fig. 9
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.
Table 1
TV comparison of the TIE algorithm and APBA. In the table, κ is the Tikhonov regularization parameter, Δ represents the sampling step-size in FT-space.

4. Discussion and Conclusion

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. (23) 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 OTFG.U.(u/M) for the geometrical unsharpness associated with the finite-size source [20]:


where Ispot([Xi w/ right arrow above]) 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]:

equation image

where operator [multiply sign in circle] denotes the convolution, A2 is the total attenuation of the imaged object, [var phi] is the spectrally averaged phase-shift caused by the object, and left angle bracket · right angle bracket 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 OTFG.U. (u/M), (ii). Replacing wavelength with the spectral-averaged left angle bracketλ2right angle bracket/left angle bracketλright angle bracket. In the same fashion, the duality transform An external file that holds a picture, illustration, etc.
Object name is nihms226357ig1.jpg defined in Eq. (6) should be modified with (i). Fourier space de-convolution of the measured intensity from OTFG.U.(u/M); (ii). A replacement of the parameter k defined in Eq.(7) with the spectral-average


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.

References and links

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.