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 ϕ by tissues is determined by

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 ϕ(r⃗) be the x-ray phase-shift caused by the imaged object, and A ₀(r⃗) the x-ray transmission, or the attenuation-map of the object. Then the Fourier transform of the x-ray intensity, $\widehat{ℱ}\left(I\right)\left(\overrightarrow{u}\right)={\int }_{{ℝ}^2}I\left(\overrightarrow{x}\right)\exp \left[2\pi i\overrightarrow{x}·\overrightarrow{u}\right]d\overrightarrow{x}$ , at position away from the object with distance R ₂, of the monochromatic point x-ray source starting at a place away from the object R ₁ distance, can be modeled by the following [9]

where I _in is the incident x-ray intensity at R ₁, λ the wavelength of the monochromatic point x-ray source and M = (R ₁ + R ₂)/R ₁ is the geometric magnification. When the FT-Space Fresnel propagator πλR ₂ u⃗ ²/M ≪ 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 πλR ₂ u⃗ ²/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 (Section3.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 A ₀ 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:

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:

with η = E _photon/m _e c ². Here we denote the photon energy of the primary x-ray beam by E _photon and m _e c ² 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 ₀ ≈ 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 ₁ from the source. We denote R ₂ as the distance from object to detector plane, M = (R ₁ +R ₂)/R ₁ the geometric magnification factor, I _in the entrance x-ray intensity at R ₁, and I_D(r⃗_D) the x-ray intensity at the detector plane. For convenience, we denote I = M ² · I_D(r⃗_D)/I _in as the normalized intensity of phase-contrast image. When the phase-attenuation duality holds, the phase map ϕ(r⃗) can be robustly retrieved from just a single projection image[14]:

where

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 ₀ as

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 ₀ = 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 ² ₀ acquired with R ₂ = 0, and the other is the acquired phase contrast image I = M ² · I_D(r⃗_D)/I _in with R ₂ > 0. Employing the acquired phase contrast image I and the duality transform Eq. (6), we will first obtain an estimate for the attenuation-component A _KN and phase map ϕ. We then rewrite Eq. (8) as

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 wavefield δAe^iϕ from R ₁ to R ₂. We can find δI = |𝔉𝔯 (δAe^iϕ)|², which is the difference between phase contrast image I and the “duality-only” counterpart I _KN = |𝔉𝔯 (A _KN e^iϕ)|². With the corrected “duality-only” image intensity $I_{\mathrm{KN}}={(\sqrt{I}+\sqrt{\delta 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 $\sqrt{I}=\sqrt{I_{\mathrm{KN}}}-\sqrt{\delta I}$ is generally valid, since it is actually a result of x-ray Fresnel diffraction and extremely smallness of hard x-ray wave-length 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 𝔉𝔯 acting on the initial wavefield, therefore the two resulted wavefields 𝔉𝔯 (A _KN exp[iϕ]) and 𝔉𝔯 (δA exp[iϕ]) would have the same resultant phases, so the resultant amplitude from the two-wave superposition is given as $\sqrt{I}=\sqrt{I_{\mathrm{KN}}}-\sqrt{\delta I}$ .

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

 

Fig. 1. Flow chart of APBA

Download Full Size | PDF

Algorithm. In an imaging experiment, assume we have performed two imaging measurements. One image (radiograph) is the attenuation image A ² ₀ acquired at SID= R ₁ , another is a normalized phase-contrast image I = M ² · I_D(r⃗_D)/I _in acquired at SID= R ₁ + R ₂. With A ² ₀ and I as well as the initial δI, usually 0, as known inputs, we first assume $I_{\mathrm{KN}}={(\sqrt{I}+\sqrt{\delta I})}^2$ . Then

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+1 − δI _k∥ is less than a predefined threshold value β, or the iteration step exceeds a predefined maximum number of iteration steps, then ϕ 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

 

Fig. 2. Example “Lena” images used in measuring the closeness between images. (a) ϕ ₁, (b) ϕ ₂, (c) ϕ ₃

Download Full Size | PDF

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 μ = ∫_Ω(g − f)d 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 ₁ norm of the gradient of the image function, and not the L ₂ norm. For two image functions f and g, we define the TV norm of g − f as the closeness of the two images:

where ∇ is the gradient operator. Since std is the form of L ₂ 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 ϕ ₁ − ϕ ₂ and ϕ ₁ − ϕ ₃ are std(ϕ ₁,ϕ ₂) = 0.206, std(ϕ ₁, ϕ ₃) = 0.472 respectively. But the image ϕ ₃ is much more closer to ϕ ₁ than ϕ ₂ 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 ϕ ₁ − ϕ ₂ and ϕ ₁ − ϕ₃ are TV(ϕ ₁, ϕ ₂) = 0.0247 and TV(ϕ ₁,ϕ ₃) = 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 ₀ 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, R ₁, and object to detector, R ₂, are set to 26 inches (0.66 m) respectively. In this way the magnification factor M = (R ₁ + R ₂)/R ₁ is equal to 2. For convenience, the incident x-ray intensity I _in at R ₁ is set to M ² (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.

 

Fig. 3. Image manifest of (a) A ² _pe,coh, (b) A ² _KN and (c) A ² ₀ when x-ray energy equals 35.5 keV.

Download Full Size | PDF

 

Fig. 4. Profiles of A ² _pe,coh, the solid lines, and A ² _KN, the dotted lines, when x-ray energy equals (a) 18.5 keV, (b) 35.5 keV and (c) 59.5 keV.

Download Full Size | PDF

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 ² _pe,coh and A ² _KN 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.

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

3.2. Comparison with the GS Algorithm

 

Fig. 5. Image representation of the inputs generated from the simulation model and Fresnel propagation. (a) the phase map ϕ; (b) the attenuation map A ² ₀; and (c) the normalized Fresnel propagated phase contrast image I with object to detector distance R ₂ = 26 in (0.66 m).

Download Full Size | PDF

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 ₂ ≈ 1 m, and converges faster for larger R ₂, 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 ₁ = R ₂ = 26 inches (0.66 m) respectively. For convenience, the incident x-ray intensity I _in at R ₁ is set to M ², where M = (R ₁ + R ₂)/R ₁ is the magnification factor. The phase map ϕ and the attenuation A ² ₀ 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 ² ₀ and the phase-contrast image I.

 

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 with the GS algorithm after 100 iterations; (c) recovered phase map with APBA after 100 iterations.

Download Full Size | PDF

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.

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 ϕ is given by

for the given normalized phase-contrast image I and attenuation image A ² ₀. Here ∇ is the gradient operator, and ∇⁻² 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.

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 ∇⁻² 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 ∇⁻² by |u⃗|²/(|u⃗|⁴ + κ ²), 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 κ is as close to the true phase ϕ 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 ϕ 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

 

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 TV measure, 0.00215226, of APBA to achieve to the TV measure, 0.00215269, of TIE; (b) recovered phase using the TIE algorithm; (c) recovered phase using APBA after 1500 iteration steps.

Download Full Size | PDF

It is L ₂ 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].

 

Fig. 8. Comparison of the TIE algorithm and APBA with noise added. (a) True phase map ϕ; (b) attenuation map A ² ₀; (c) the normalized phase contrast image I; (d) recovered phase map with the TIE algorithm, no Tikhonov regularization is used; (e) recovered phase map with the TIE algorithm with Tikhonov regularization; (f) recovered phase map with APBA after 10 iteration steps. In the simulation, the acquired data is assumed to have a level of δ_b = 0.03% detector noise and one pixel misalignment between A ² ₀ and I horizontally.

Download Full Size | PDF

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 ² ₀ 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 σ ² = 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/N_b corresponding to the direct exposure area (where A ² = 1) outside the object in background. The Poisson statistics dictates that N_b = 1/δ ² _b 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 ${\delta }_b=\sigma P_b\sqrt{P_b}$ is the statistical coefficient of variation in absorbed photon numbers in background. ${\delta }_{\mathrm{oi}}=\sigma P_{\mathrm{oi}}\sqrt{P_{\mathrm{oi}}}$ , 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 A ² ₀ 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 A ² ₀ although they have the same background noise level.

Three simulation tests are performed in the comparison: case1: assume the acquired data A ² ₀ 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.

4. Discussion and Conclusion

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.

View Table

 

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.

Download Full Size | PDF

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 theWigner 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.(u⃗/M) for the geometrical unsharpness associated with the finite-size source [20]: ${\mathrm{OTF}}_{G.U.}\left(\frac{\overrightarrow{u}}{M}\right)=\frac{\int I_{\mathrm{spot}}\left(\overrightarrow{\xi }\right)\exp [i2\pi \overrightarrow{\xi }·\frac{\left(M-1\right)\overrightarrow{\mathbf{u}}}{M}]d\overrightarrow{\xi }}{\int I_{\mathrm{spot}}\left(\overrightarrow{\xi }\right)d\left(\overrightarrow{\xi }\right)}$

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]: $I(\overrightarrow{\mathbf{r}};R_1+R_2)=\frac{I_{\mathrm{in}}}{M^2}{\widehat{ℱ}}^{-1}\left({\mathrm{OTF}}_{G.U.}\left(\frac{\overrightarrow{\mathbf{u}}}{M}\right)\right)\otimes \{A^2\left(\frac{\overrightarrow{\mathbf{r}}}{M}\right)-\frac{R_2〈{\lambda }^2〉}{2\pi M〈\lambda 〉}\nabla ·\left(A^2\nabla \varphi \left(\frac{\overrightarrow{\mathbf{r}}}{M}\right)\right)\},$

where operator ⊗ denotes the convolution, A ² 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.(u⃗/M), (ii). Replacing wavelength with the spectral-averaged 〈λ ²〉/〈λ〉. In the same fashion, the duality transform 𝔇 defined in Eq. (6) should be modified with (i). Fourier space deconvolution of the measured intensity from OTF_G.U.(u⃗/M); (ii). A replacement of the parameter k͂ defined in Eq.(7) with the spectral-average $〈\tilde{k}〉=\frac{r_e{R}_2}{2\pi M}·〈\frac{{\lambda }^2}{{\sigma }_{\mathrm{KN}}}〉.$ .

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.