1. Introduction

There has been a surge of interest recently in the theories and algorithms for in-line x-ray phase imaging and phase retrieval [1–12] due to their potential advantages over traditional purely attenuation-based x-ray imaging in clinical applications like mammography and general radiography. Several groups have presented [5, 8, 11, 12] their formulations based on the paraxial Fresnel-Kirchoff diffraction integral or the Transport of Intensity Equation (TIE), which is an approximation of the diffraction integral in near field [11, 13]. To facilitate algorithm implementation of direct phase retrieval, the formulations were in a “linearized” form so that the phase term can be solved explicitly. These theories have been used to retrieve the phase maps of several simple models such as (a) “weak” objects (where both attenuation and phase shift are small), (b) “pure-phase” objects (where attenuation is constant) and (c) “homogeneous” objects (where the phase and the logarithm of intensity are proportional to each other). However, studies on general objects with nontrivial and uncorrelated distributions of phase and attenuation are much less successful so far with these methods.

Wu and Liu presented a more general and concise formulation [5] and later extended it in phase-space [8] to treat polychromation and finite focal spot size of x-ray sources and finite resolution of detectors. This approach was also linearized by neglecting two terms of the formula.

In this communication, we will examine the effect of the linearization and point out that the terms could be kept for a more quantitative and more physically reasonable solution. Based on this result, we present an iterative algorithm to solve the nonlinear problem. Simulation studies show that, when dealing with complex objects, this iterative algorithm can give more quantitative result than the linearized one. The algorithm also has high computational efficiency and scalability.

2. Theoretical background

A “thin sample” [5] in x-ray imaging can be described by its optical transmission function T(x,y) = A(x,y)exp(iϕ(x,y)). The objective of phase retrieval is to obtain the phase mapϕ(x,y) from an attenuation-based image I ₁ taken with M = 1 and a phase-contrast image I ₂ taken with M > 1. Here M = (R ₁ + R ₂)/R ₁ is the geometrical magnification, R ₁ and R ₂ are the source-object-distance (SOD) and the object-image-distance (OID).

For simplicity, we use one-dimensional notions of formulas in this communication. Using the “back-projected image” I _bp(x) = M ² I(Mx), the formulation of phase imaging by Wu and Liu [5] can be written as:

where I ₀ is the incident x-ray intensity, λ is the wavelength of the x-ray, u is the spatial frequency, $\widehat{ℱ}$ [∙] means Fourier transform. This formulation is obtained through the moderate variation conditions of the attenuation and phase maps:

which are of the forms of Taylor expansions and can easily be satisfied in clinical imaging due to the extremely short wavelength λ of the x-ray and the limited spatial resolution u of the detectors. As has been estimated for a typical mammography system [5], λR ₂ u/M is less than 0.62 micron, which is much smaller than the typical size of the breast structures to be imaged. In such cases, the moderate variation conditions can be met. This formula Eq. (1) considers the most general cases in clinical applications in a very concise form, providing high flexibility in the implementation of phase retrieval algorithms. The formulation was further extended to include the effects of polychromatic x-ray source with finite focal spot size and real detectors with finite spatial resolution [14] by multiplying the right-hand side of Eq. (1) with the OTFs of the x-ray source, OTF_G.U, and the detector, OTF_det.

Eq. (1) is a nonlinear formula and is difficult for analytical solution. To solve ϕ(x,y) from Eq. (1), it was linearized regarding the Fourier transform of A ²(x,y)ϕ(x,y) as [5]

from which $\widehat{ℱ}$[A ²(η)ϕ(η)] and hence ϕ can be solved and calculated explicitly.

3. The iterative algorithm

Now we will examine the linearization of Eq. (1) by comparing the order of magnitude of each term. For simplicity, we denote the four terms as T ₁ through T ₄. The assumption of A(η) ≈ A(η ± λR ₂ u/M) used for the linearization [5] could be valid, however, it cannot assure that T ₂ ≫ T ₄ since sin (λR ₂ u ²/M) can be much smaller than cos (λR ₂ u ²/M). T ₁ through T ₄ are dependent on the actual distribution of both A(η) and ϕ(η) and therefore no general conclusion can be made about their magnitudes for all cases. Here we just give some very rough estimations.

One can see that, since λR ₂ u ²/M ≫ 1 holds true for typical clinical applications, T ₁ through T ₄ are approximately proportional to λ⁰ λ¹, λ² and λ¹, respectively, which leads to the estimation that

Furthermore, it is known that adding a constant to the phase map will not change the physics and hence the obtained image, which means

From this result, one can also estimate that T ₂ and T ₄ are of the same order, at least for some cases. Therefore, in phase retrieval studies, of the four terms of Eq. (1), T ₃ can be neglected while T ₄ should be kept. With the effects of the x-ray focal spot size and the detector resolution taken into account as presented by Wu and Liu [8], the formula becomes

The discussion above also leads to the conclusion that, besides numerical exactness, keeping T ₄ is a requirement of physical correctness. Without T ₄, the image will be dependent on the choice of the zero point of ϕ(η), which is physically incorrect, while keeping T ₄ in the formula can avoid the ambiguity.

If T ₄ is kept, the formula remains nonlinear and is not explicitly solvable, and therefore, direct retrieval is no longer applicable. To treat this problem, we present an iterative retrieval algorithm as follows:

Then repeat Step 3–4 until a steady distribution of A ² ϕ is obtained. Dividing A ² ϕ by A ² gives the phase map ϕ

One problem in the algorithm is how to determine the zero-frequency amplitude of the retrieved A ² in Step 4, which cannot be obtained by dividing 2 sin (πλR ₂ u ₂/M) since sin(0) = 0. This problem was solved by choosing the zero point of phase map ϕ and hence that of A ² ϕ. Practically, one can choose a “region of zero phase” (RZP) in the object plane, in which both ϕ and A ² ϕ are assumed to be zero. The choice of RZP is arbitrary as long as ϕ is approximately a constant (not necessarily zero) in it. In the algorithm, after each Step 4, a constant equal to the average value in RZP is subtracted from A ² ϕ Thus, the zero-frequency amplitude is determined. It should be pointed that, though the choice of RZP does not affect the distribution of ϕ, it does on A ² ϕ because A ²(ϕ + c) will generally differ from A ² ϕ for a constant c.

It is worth noting that, in most x-ray images, there is “open field”—an area not occupied by the object being imaged—within which the phase map is relatively uniform. The open field can favorably serve as the RZP. And in this case, the retrieved phase map exactly equals the phase shift cased by the object itsef.

4. Tests and results

The above iterative algorithm was tested using simulated images of a “virtual sample”. The simulation is implemented by calculating the paraxial Fresnel-Kirchhoff diffraction integral directly. For simplicity, an ideal monochromatic point x-ray source and an ideal detector were assumed in the simulations. To test the quantitativeness of the algorithm, we used a sample more complicated than those simplified ones studied before, with distinct and nontrivial distributions of both attenuation and phase maps.

The parameters for the simulation are: kVp= 20keV, pixel pitch= 50μm, R ₁ = R ₂ = 1m. The virtual sample is of an attenuation map “cameraman” as shown in Fig. 1(a) and a phase map “Lena” as shown in Fig. 1(b). With the distinct attenuation and phase map, one can see more evidently the effectiveness of the retrieval algorithm. The RZP is chosen as the 100 by 100 square on the upper-left corner of the sample. In the RZP, ϕ = 5.

 

Fig. 1. The attenuation and phase maps of the virtual sample for the simulation. The size of both maps is 1024 by 1024. The attenuation map is of average value 0.6455 and standard deviation 0.2257. The phase map is of average value 9.3203 and standard deviation 4.8829.The RZP is chosen as the 100 by 100 square on the upper-left corner of the sample. In the RZP, ϕ = 5.

Download Full Size | PDF

The attenuation-based map I ₁ and the phase-contrast map I ₂ obtained by the simulation program are shown in Fig. 2. One can see that, due to the short wavelength of the x-ray and hence weak diffraction, the phase contrast effects are almost unnoticeable in I ₂, but can be seen in the difference image M ² I ²-I ₁. The standard deviation (STD) of the difference image is 2.8×10^-4, three orders of magnitude smaller than that of the attenuation map.

 

Fig. 2. The simulated attenuation-based and phase-contrast images of the virtual sample. Due to the short wavelength of the x-ray and hence weak diffraction, the phase contrast effects are almost unnoticeable in I ₂, but can be seen in the difference image M ² I ₂ -I ₁.

Download Full Size | PDF

For comparison, we first did the retrieval using the linearized algorithm. The results are shown in Fig. 3. One can see that the retrieved map of A ² ϕ contains information of both attenuation and phase, however, it is clearly not the product of the two. Moreover, as discussed above, there exists ambiguity in the retrieval and the phase map cannot be uniquely determined. Direct division of the obtained A ² ϕ by A ² generally tends to give a phase map with bad or even wrong distribution as shown in Fig. 3(b).

 

Fig. 3. The retrieved phase map using the linearized formula. The retrieved map of A ² ϕ is clearly not the product of the two. And due to the ambiguity in the retrieval, the phase map cannot be uniquely determined. Direct division of (a) gives a phase map (b).

Download Full Size | PDF

The results of the iterative algorithm after 20 iterations are shown in Fig. 4. The retrieved phase map is very good except for some remnants of attenuation map where sharp attenuation contrast exists. This is partly due to the deficiency of the algorithm currently used to calculate the gradient of the attenuation map and partly due to the violation of the assumption Eq. (3) and may be improved in further studies as will be discussed later.

 

Fig. 4. The retrieved phase map using the iterative algorithm. The small remnants of attenuation on ϕ is due to the deficiency of the algorithm for calculating the gradient of the attenuation map and the violation of the assumption Eq. (3).

Download Full Size | PDF

To test the convergence of the algorithm, we used the the difference image Δϕ between the retrieved map and the original map used for simulation. The STD of Δϕ is plotted against the iteration count in Fig. 5.

 

Fig. 5. The standard deviation of the difference map versus the iteration count. The STD increases at first because of the improperly chosen initial distribution A2f = 0, but converges quickly and levels off after iteration 13.

Download Full Size | PDF

The arbitrarily chosen initial distribution A ² ϕ = 0 leads to the increase of the STD for the first 4 iterations. Even though, the algorithm is robust enough to find a convergent solution quickly after 13 iterations. The algorithm is also of high efficiency. Actually, the 20 iterations was done on an ordinary PC in several minutes with large scope for further optimization.

For more quantitative comparison, we also did one dimensional tests using simulations. The attenuation and phase maps of the one dimensional “virtual object” are taken from the 500th row of the images shown in Fig. 1(a) and (b). The RZP is the first 10 pixels from left. Other parameters are the same with the simulations above. One-dimensional formulas are used to obtain the result, which is shown in Fig. 6.

 

Fig. 6. One dimensional tests of the iterative retrieval algorithm. The attenuation and phase maps of the “virtual object” are taken from the 500th row of the images shown in Fig. 1(a) and (b), respectively.

Download Full Size | PDF

It is evident that the result by the iterative algorithm is more quantiative as compared to that by the linear algorithm. Moreover, there are spikes in the curve obtained by the linear algorithm at places the attenuation changes sharply, whereas the iterative algorithm is much less affected.

5. Discussions and conclusions

One can see that, with the inital distribution A ² ϕ = 0, the linearized formula is actually the same as the first iteration of the iterative approach. Furthermore, for “pure-phase” objects with constant attenuation, the iterative formula degrades to the linear one. Therefore, the iterative algorithm can be seen as an extension to the linearized algorithm, with higher precision and flexibility. Nevertheless, from another point of view, it is more than just an extension. Traditional iterative phase retrieval approaches used in crystallography such as the Gerchberg-Saxton (GS) scheme [15] are proved to be less deterministic, robust and computationally efficient, whereas popular linear retrieval algorithms based on the paraxial Fresnel-Kirchoff diffraction integral are not so successful in handling complex objects. Based on an extensive and general theoretical formulation for in-line phase imaging, our iterative algorithm has the merits of both quantitativeness and efficiency, therefore has great potential for complicated samples in clinical applications.

The iterative algorithm also has high scalability to treat even more complex conditions than has been discussed, i.e., more rapidly varying attenuation maps not conforming to Eq. (3). During the derivation of Eq. (1), only the leading terms of the Taylor expansions were kept in Eq. (3) to make possible the linearization. With the iterative nonlinear approach, however,one can further enhance the acuracy of the algrothm by keeping more terms to better describe the attenuation and phase maps. In this way, we can hopefully reduce the remnants of the attenuation map in Fig. 4(b). This will be a topic of further studies.

Though the algorithm proves robust to the initial value of the iteration, problems remain regarding its robustness agaist the system noise. In the simulated system, the sum of the phase-related terms is approximately 3 orders smaller than that of the attenuation-based map. Generally, the level of quantum noise in practical systems is larger than that. Special considerations need be taken for both the design of the whole imaging system and the retrieval algorithm to handle the interference of the system noise. This is another important topic under investigation and will be reported in a future publication.