Noise propagation in x-ray phase-contrast imaging and computed tomography

All experimentally acquired images contain noise. The nature of the noise can be different and strongly depends on the properties of the radiation source and the detection system. A detailed review of this topic can be found in [9]. In the present work, we ignore the issue of detector efficiency as a function of x-ray energy. We also restrict our consideration to noise associated with Poisson photon counting statistics, thus implicitly assuming a Poisson x-ray source and a photon counting detector in the imaging system.

Also, we assume that noise in projection images of an object is approximately spatially stationary, i.e. the statistical properties of noise are spatially shift-invariant. Although in reality the noise in any non-uniform and/or finite image is never strictly stationary, it can be very accurately treated as stationary as long as the contrast in the image is weak and the size of the image is significantly bigger than the width of the noise autocorrelation function. If however the contrast in the image is strong, but the intensity, as the function of the spatial location in the image, is slowly varying on the length scale of the noise autocorrelation function, then the corresponding noise can be treated as quasi-stationary [9].

An important characteristic of random noise is NPS that, in the case of one-dimensional (1D) noise n(x), is defined as follows [9, section 8.2.5],

$$\begin{equation} W_{n} (u)\stackrel{\Delta}{=} \lim\limits_{L\to \infty} \frac{1}{L}\left\langle {\left| {\int_{-L/2}^{L/2} {{\rm d}xn(x)\exp (2\pi {\rm i}ux)}} \right|^{2}} \right\rangle , \end{equation} \tag{ 2.1 }$$

where the triangular brackets denote an ensemble average. It should be noted that W_n(u) is defined for non-stationary as well as stationary random noise. However, some transformations of NPS are more straightforward in the case of stationary noise. For instance, convolution of a 1D image I(x) containing noise with some point-spread function, which can be conveniently presented as filtering in Fourier space, $\hat{{I}}_{T} (u)=\hat{{I}}(u)\,\hat{{T}}(u)$, results in the following simple transform of the corresponding NPS [9]:

$$\begin{equation} W_{n,T} (u)=W_{n} (u)\,\vert \hat{{T}}(u)\vert^{2}. \end{equation} \tag{ 2.2 }$$

Here and throughout the paper, the triangular cap ^∧ designates the Fourier transform.

In the case of quasi-stationary noise in a non-uniform image, its NPS (in 1D case) depends on the spatial coordinate x via a scaling factor, W_n(u; x) = I(x)W_n,0(u), where W_n,0(u) is the NPS in a sample-free (uniform) image and I(x) is the image intensity (divided by the corresponding flat field).

Two other important characteristics of random noise are its mean E_n and variance ${\rm std}_{n}^{2}$, that, in the 1D case, can be defined as follows:

$$\begin{equation} E_{n} \equiv \lim\limits_{L\to \infty} \frac{1}{L}\left\langle {\int_{-L/2}^{L/2} {{\rm d}xn(x)}} \right\rangle \end{equation} \tag{ 2.3 }$$

and

$$\begin{equation} {\rm std}_{n}^{2} \equiv \lim\limits_{L\to \infty} \frac{1}{L}\left\langle {\int_{-L/2}^{L/2} {{\rm d}x[n(x)-E_{n} ]^{2}}} \right\rangle . \end{equation} \tag{ 2.4 }$$

For any continuous 1D noise distribution with zero mean (we will implicitly assume this hereafter), its variance ${\rm std}_{n}^{2}$ is related to the corresponding NPS W_n(u) as follows¹:

$$\begin{equation} {\rm std}_{n}^{2} =\int_{-\infty}^{+\infty} {{\rm d}uW_{n} (u)} . \end{equation} \tag{ 2.5 }$$

Taking into account that all real images are finite and discrete, their NPS is discrete and periodic, in the Fourier space, with the period 2U_m = h⁻¹ defined by the image pixel size h. Assuming that the image size (in pixels) is large enough, so that the NPS is finely sampled, one can estimate the variance of noise using the following integral (which is more convenient for analytical analysis presented below than the summation over discrete Fourier spectrum),

$$\begin{equation} {\rm std}_{n}^{2} \cong \int_{-U_{m}}^{+U_{m}} {{\rm d}uW_{n} (u)} . \end{equation} \tag{ 2.6 }$$

In particular, for uncorrelated (white) noise, equation (2.6) gives the following simple result:

$$\begin{equation} W_{n,{\rm wht}} =h\,{\rm std}_{n,{\rm wht}}^{2} . \end{equation} \tag{ 2.7 }$$

Similarly to equation (2.6), for any discrete two-dimensional (2D) noise distribution with zero mean, the variance of random noise is related to the corresponding NPS W_n(u, v) as follows:

$$\begin{equation} {\rm std}_{n}^{2} \cong \int_{-U_{m}}^{+U_{m}} {{\rm d}u} \int_{-U_{m}}^{+U_{m}} {{\rm d}vW_{n} (u,v)} . \end{equation} \tag{ 2.8 }$$

PBI [10, 11] is the simplest PCI method that does not require any optical elements between x-ray source, object and detector. In this paper we restrict our consideration to the case of a plane monochromatic wave incident on the object. However, the results obtained here can be easily extended to the case of a divergent polychromatic beam (e.g. spherical wave) incident on the object. The former case is usually attributed to synchrotron x-ray sources while the latter case is typical for (microfocus) laboratory x-ray sources.

In this section we review two widely used phase-retrieval approaches, based on the transport-of-intensity equation (TIE) [12] and the CTF [13]. We apply these two approaches to monomorphous objects, i.e. the objects in which the ratio γ ≡ δ/β is spatially invariant. This property results in the following linear dependence between the phase shift φ and the attenuation function B induced by the object (the latter is related to object's contact image I₀, I₀ = exp(−2B)):

$$\begin{equation} \varphi(x,y) =-\gamma B(x,y). \end{equation} \tag{ 2.9 }$$

Hereafter, we implicitly assume that all projection images of an object (including contact and phase-contrast images) are divided by the corresponding flat fields (i.e. images collected without the object), or, more precisely, we assume a plane incident wave with unit intensity.

2.2.1. TIE-HOM.

TIE approach [12], applied to a monomorphous object (see equation (2.9)), relates a contact image I₀ of the object and a phase-contrast image I_z collected at some distance z from the object, downstream in the propagation direction of the x-ray beam [14]:

$$\begin{equation} \hat{{I}}_{z} (u,v)=\left[ {1+\pi \gamma \lambda z(u^{2}+v^{2})} \right]\;\hat{{I}}_{0} (u,v), \end{equation} \tag{ 2.10 }$$

where λ is the x-ray wavelength. Following [14], we will call this equation TIE-HOM.

Equation (2.10) allows one to reconstruct the contact image I₀ of the object using only one phase-contrast image I_z. Using equation (2.10) and a 2D analogue of equation (2.2), NPS in the reconstructed contact image, W_n,0, is related to NPS W_{n, z} in the corresponding phase-contrast image as follows:

$$\begin{equation} W_{n,0} (u,v)=\frac{W_{n,z} (u,v)}{[1+\pi \gamma \lambda z(u^{2}+v^{2})]^{2}}. \end{equation} \tag{ 2.11 }$$

In turn, NPS in the phase-contrast image can be written as follows:

$$\begin{equation} W_{n,z} (u,v)={\rm MTF}_{{\rm det}}^{2} (u,v)W_{n,{\rm wht}} , \end{equation} \tag{ 2.12 }$$

where MTF_det is the modulation transfer function (MTF) of the detector and W_n,wht is the NPS of uncorrelated/white noise, which in view of equation (2.8) can be expressed in terms of its standard deviation, std_n,wht, as follows:

$$\begin{equation} W_{n,{\rm wht}} =h^{2}{\rm std}_{n,{\rm wht}}^{2} =h^{2}/N_{{\rm ph,pxl}} . \end{equation} \tag{ 2.13 }$$

Hereafter we assume that the only source of noise in projection images is the photon counting noise satisfying Poisson statistics. Hence, ${\rm std}_{n,{\rm wht}}^{2} =1/N_{{\rm ph,pxl}}$, where N_ph,pxl is the number of photons per detector pixel (recall that all projection images are assumed to be divided by the corresponding flat field; and the phase contrast is assumed to be weak).

2.2.2. CTF-HOM.

CTF approach has been originally developed by Guigay [13] for weakly absorbing objects (so that |B| ≪ 1) with the phase φ satisfying the following condition, |φ(x, y) − φ(x − λzu, y − λzv)| ≪ 1, which can be interpreted that the phase is either slowly varying or small. Under these assumptions, the Fourier transform of the intensity I_z in a phase-contrast image can be written as follows [13]:

$$\begin{eqnarray} &&\fl\hat{{I}}_{z} (u,v)\cong \delta_{{\rm D}} (u)\delta_{{\rm D}} (v)-2\hat{{B}}(u,v)\cos [\pi \lambda z(u^{2}+v^{2})] \nonumber\\ &&\quad+\,2\hat{{\varphi}}(u,v)\,\sin [\pi \lambda z(u^{2}+v^{2})], \end{eqnarray} \tag{ 2.14 }$$

with δ_D(x) denoting the Dirac delta function.

For a monomorphous object, for which equation (2.9) is satisfied, equation (2.14) simplifies,

$$\begin{eqnarray} &&\fl\hat{{I}}_{z} (u,v)\cong \hat{{I}}_{0} (u,v)\sqrt {\gamma^{2}+1} \,\sin [\pi \lambda z(u^{2}+v^{2})+\arctan (\gamma^{-1})].\nonumber\\ \end{eqnarray} \tag{ 2.15 }$$

We will call this equation CTF-HOM. In order to avoid zeros of the sine function in equation (2.15), the reconstruction of the contact intensity I₀ is usually carried out using the following regularized formula:

$$\begin{eqnarray} &&\fl\hat{{I}}_{0} (u,v)\nonumber\\ &&\fl\quad\approx \frac{\hat{{I}}_{z} (u,v)\,\sin [\pi \lambda z(u^{2}+v^{2})+\arctan (\gamma^{-1})]}{\,\sqrt {\gamma^{2}+1} \,\{\sin ^{2}[\pi \lambda z(u^{2}+v^{2})+\arctan (\gamma^{-1})]+\,\varepsilon (u,v)\}},\nonumber\\ \end{eqnarray} \tag{ 2.16 }$$

where, in the case of positive γ, the regularization parameter ε is

$$\begin{eqnarray} \fl\varepsilon (u,v)=\left\{{\begin{array}{@{}ll@{}} 0, & \pi \lambda z(u^{2}+v^{2})+\arctan (\gamma^{-1})<\pi /2, \\ \varepsilon_{0} ,&{\rm otherwise.} \end{array}} \right.\nonumber\\ \end{eqnarray} \tag{ 2.17 }$$

Here ε₀ ≪ 1 is a small regularization constant.

Accordingly, NPS in the reconstructed contact image is written as follows:

$$\begin{eqnarray} &&\fl W_{n,0} (u,v)\nonumber\\ &&\fl\quad\cong \!\frac{W_{n,z} (u,v)\,\sin^{2}[\pi \lambda z(u^{2}+v^{2})+\arctan (\gamma^{-1})]}{\,(\gamma^{2}+1)\{\,\sin^{2}[\pi \lambda z(u^{2}+v^{2})+\arctan (\gamma^{-1})]+\,\varepsilon (u,v)\}^{2}}.\nonumber\\ \end{eqnarray} \tag{ 2.18 }$$

Equations (2.12) and (2.13) can be used to estimate NPS in a phase-contrast image, W_{n, z}.

Although a rigorous analysis of the validity of equation (2.15) is beyond the scope of this paper, nevertheless, some indications of this are given below. First of all, it should be noted that TIE is valid for any slowly varying wave field in the object plane (with no additional limitations on the attenuation magnitude) such that πλz(u² + v²) ≪ 1 (see, for instance, [15]). As a result, equation (2.10) can be formally transformed into equation (2.15) using the following substitutions: sin[πλz(u² + v²)] ≈ πλz(u² + v²) and cos[πλz(u² + v²)] ≈ 1.

For comparison, CTF approach (in its original formulation by Guigay [13]) imposes more stringent requirement on the attenuation magnitude, |B| ≪ 1, but less stringent requirement on the phase, in the form of Guigay's condition above. Alternatively, in [15] it has been shown that equation (2.15) also describes phase-contrast images of a monomorphous object whose complex transmission function q(x, y) can be presented as follows, $q(x,y)=\bar{{q}}(x,y)\,[1+\psi (x,y)]$. Here $\bar{{q}}(x,y)$ is a slowly varying component satisfying the TIE validity condition, namely, that it can be linearized in the vicinity (with the linear size of λ zu) of any point (x, y) in the object plane; and ψ(x, y) is a small distortion, |ψ(x, y)|≪ 1, so that the component 1 + ψ(x, y) satisfies the so-called weak object approximation [13]. This composite approach is called TIE+Born in [15].

In this section we review a widely used PCI method that utilizes a perfect analyser crystal, the so-called analyser-based imaging (ABI) [16–18]. We restrict our consideration to a phase-retrieval algorithm which is based on the geometrical optics approximation (GOA) [19]. According to this approximation, intensity distribution I(x, y) in the phase-contrast image is related to the phase φ and the attenuation function B of the object as follows:

$$\begin{equation} I_{\omega} (x,y)=I_{0} \,R(\omega -k^{-1}\partial_{x} \varphi ), \end{equation} \tag{ 2.19 }$$

where I₀ ≡ exp(−2B) is the contact image, R(ω) is the rocking curve of the analyser crystal and ω is the deviation of the analyser crystal from the exact Bragg angular position, k = 2π/λ. In equation (2.19) we implicitly assume (through the use of the partial derivative ∂_xφ) that the diffraction plane of the analyser crystal is perpendicular to y-axis and parallel to x-axis.

Let us make the following assumptions [20]:

• (1)  
  the object is monomorphous, equation (2.9);
• (2)  
  the rocking curve of the analyser crystal is well approximated by a Gaussian function,

  $$\begin{equation} R(\omega )=R_{0} \exp (-\omega^{2}/\sigma^{2}); \end{equation} \tag{ 2.20 }$$
• (3)  
  local angular deviations, k⁻¹∂_xφ, of the wavefront incident on the analyser crystal are small compared to its deviation angle ω,

  $$\begin{equation} \vert k^{-1}\partial_{x} \varphi \vert \ll \vert \omega \vert . \end{equation} \tag{ 2.21 }$$

Substituting equation (2.20) into equation (2.19) and neglecting, in accordance with equation (2.21), the quadratic term, one obtains

$$\begin{equation*} \ln I_{\omega} \cong \ln [I_{0} R(\omega )]+2\sigma^{-2}\omega \,k^{-1}\partial_{x} \varphi . \end{equation*}$$

The second term on the rhs of the last equation is the phase contrast which is usually small compared to one (especially when equation (2.21) is satisfied and the deviation ω of the analyser crystal is not large compared to the rocking curve width), so that one can rewrite this term as follows: 2σ⁻²ω k⁻¹∂_xφ ≅ ln[1 + 2σ⁻²ω k⁻¹∂_xφ], and hence

$$\begin{equation*} I_{\omega} \cong I_{0} R(\omega )\,[1+2\sigma^{-2}\omega \,k^{-1}\partial _{x} \varphi ]. \end{equation*}$$

Using equation (2.9) and noting that −2I₀∂_xB = ∂_xI₀, one then obtains

$$\begin{equation} I_{\omega} \cong R(\omega )\{I_{0} \,+\,\gamma \sigma^{-2}\omega k^{-1}\partial_{x} I_{0} \}. \end{equation} \tag{ 2.22 }$$

Taking the Fourier transform of equation (2.22), one then has

$$\begin{equation} \hat{{I}}_{0} (u,v)\cong \frac{\hat{{I}}_{\omega} (u,v)}{R(\omega )\{1-2\pi {\rm i}\gamma \sigma^{-2}\omega k^{-1}u\}}. \end{equation} \tag{ 2.23 }$$

It should be noted that the same form of solution for I₀ as in equation (2.23) could be obtained if, instead of assuming Gaussian shape for the rocking curve of the analyser crystal, one would expand the rocking curve function into the Taylor series,

$$\begin{eqnarray*}&&\fl R(\omega -k^{-1}\partial_{x} \varphi )\approx R(\omega )-k^{-1}\partial_{x} \varphi \,{R}'(\omega )\nonumber\\ &&+(1/2)(k^{-1}\partial_{x} \varphi )^{2}{R}''(\bar{{\omega}}), \end{eqnarray*}$$

and then neglect the quadratic term [21]. Then, applying equation (2.9), one would have the following expression for the reconstructed contact image,

$$\begin{eqnarray} &&\fl\hat{{I}}_{0} (u,v)\cong \frac{\hat{{I}}_{\omega} (u,v)}{R(\omega )\{1+\pi {\rm i}\,\gamma [{R}'(\omega )/R(\omega )]\,k^{-1}u\}}. \end{eqnarray} \tag{ 2.24 }$$

Note that equation (2.24) is applicable to a linear region of any rocking curve. In particular, assuming a Gaussian shape of the rocking curve in this equation, one immediately obtains equation (2.23) since for the Gaussian shape of R, as in equation (2.20), R'(ω)/R(ω) = −2σ⁻²ω.

For our further analysis, it is convenient to introduce in equation (2.23) a dimensionless parameter c,

$$\begin{equation} c=2\pi \,\gamma (\omega /\sigma )\,(\Delta \omega /\sigma ), \end{equation} \tag{ 2.25 }$$

where Δω ≡ k⁻¹U_m = λ/(4πh) corresponds to the maximum measured diffraction angle. Equations (2.23), (2.25) and (2.2) allow one to write NPS for the reconstructed contact image,

$$\begin{equation} W_{n,0} (u,v)\cong \frac{W_{n,\omega} (u,v)}{R^{2}(\omega )[1+c^{2}\,(u/U_{\max} )^{2}]}. \end{equation} \tag{ 2.26 }$$

It should be mentioned that for a fixed incident photon fluence (and, therefore, for a fixed radiation dose delivered to the object), the photon counting noise in a phase-contrast image depends on the analyser-crystal deviation angle, namely, W_n,ω ∼ R(ω). Therefore, under the assumption on a constant photon fluence, the NPS in the reconstructed contact image, in view of equation (2.26), is inversely proportional to the rocking curve value, W_n,0 ∼ R⁻¹(ω).

Recall that in this paper we restrict our analysis to the case of a plane monochromatic wave incident on the object. We also assume that the projection approximation can be applied to calculate the attenuation and the phase shift acquired by the wave upon transmission through the object, i.e.

$$\begin{eqnarray} &&\fl B_{\theta} ({x}',y)=k(\mbox{\bit{P}}_{\theta} \beta )({x}',y), \quad \varphi_{\theta} ({x}',y)=-k(\mbox{\bit{P}}_{\theta} \delta )({x}',y),\nonumber\\ \end{eqnarray} \tag{ 2.27 }$$

where the x-ray transform (projection) [22, p 9] of a real-valued function f(x, y, z) is defined by $(\mbox{\bit{P}}_{\theta} f)({x}',y)=\int_{-\infty}^{+\infty} {\int_{-\infty}^{+\infty} {{\rm d}x\,{\rm d}zf(x,y,z)\delta_{{\rm D}} ({x}'-x\sin \theta -z\cos \theta )}}$, and rotation is carried out around the y-axis. Given projections (P_θf)(x',y) for all view angles θ ∈ [0,π ], the function f(x, y, z) can be uniquely reconstructed using the FBP reconstruction formula [22, p 102],

$$\begin{eqnarray} &&\fl f(x,y,z)=\int_0^\pi \int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} {\rm d}\theta \,{\rm d}\xi \,{\rm d}\eta {\rm e}^{-2\pi {\rm i}\,[\xi (x\sin \theta +z\cos \theta )+\eta y]} \nonumber\\ &&\times \vert \xi \vert (P_{\theta} f)^{\wedge}(\xi ,\eta ) .\nonumber\\ \end{eqnarray} \tag{ 2.28 }$$

If CT reconstruction of the linear attenuation coefficient μ = 2kβ is carried out using FBP reconstruction formula, equation (2.28), from a finite (but sufficiently large) number of continuous projections then 2D NPS W_n,μ(u, v) in a reconstructed continuous slice is related to 1D power spectrum W_n,0(ξ) of stationary noise in an individual contact projection as follows [23, 24]:

$$\begin{equation} W_{n,\mu} (u,v)=W_{n,\mu} (\xi )=\frac{\pi}{N_{{\rm a}}}\xi W_{n,0} (\xi ), \end{equation} \tag{ 2.29 }$$

where ξ is the radial spatial frequency defined as $\xi =\sqrt {u^{2}+v^{2}}$, N_a is the number of equidistant angular positions over the angular interval of 180° at which the CT projection data are acquired. It is worth mentioning the different dimensionality of W_n,μ(L⁰) and W_n,0(L¹). The variance of noise in the reconstructed (continuous) slice then can be calculated using a 2D analogue of equation (2.8):

$$\begin{equation} {\rm std}_{n,\mu}^{2} =2\pi \int_0^\infty {{\rm d}\xi \xi W_{n,\mu} (\xi )} . \end{equation} \tag{ 2.30 }$$

When deriving an expression for the variance of noise in a CT slice reconstructed from discrete projections one can use one of the following two alternative approaches. In the first approach proposed in [25] analysis is carried out in Fourier space (continuous). In this case, one should explicitly take into account the periodic structure of NPS of discrete projections and the filtering effect of interpolation during the backprojection step. In the second approach used e.g. in [26, 27] one directly operates with discrete Fourier transforms. Using either of the two approaches, one can easily obtain the following general expression for the variance of noise in a CT slice:

$$\begin{equation} {\rm std}_{n,\mu}^{2} =2\pi \int_0^{U_{m}} {{\rm d}\xi \,\xi W_{n,\mu} (\xi )} , \end{equation} \tag{ 2.31 }$$

where an explicit form of the NPS depends on the used Fourier filter G(ξ) (e.g. ramp, Hanning etc) [24] and on the interpolation scheme utilized in the backprojection step. For example, in the case of nearest-neighbour (NN) interpolation and linear interpolation (LI), the corresponding expressions for the NPS are written as follows:

$$\begin{equation} W_{n,\mu}^{{\rm NN}} (\xi )=\frac{\pi}{N_{{\rm a}}}\xi^{-1}G^{2}(\xi )W_{n,0} (\xi ), \end{equation} \tag{ 2.32 }$$

$$\begin{equation} W_{n,\mu}^{{\rm LI}} (\xi )=\frac{\pi}{N_{{\rm a}}}\xi^{-1}G^{2}(\xi )\frac{2+\cos (2\pi h\xi )}{3}\,W_{n,0} (\xi ). \end{equation} \tag{ 2.33 }$$

Propagation of the variance of random noise from a phase-contrast image to the corresponding contact image obtained with the use of TIE-HOM phase retrieval is described by equation (2.8) with the NPS W_n,0 given by equations (2.11)–(2.13),

$$\begin{equation} {\rm std}_{n,0}^{2} ={\rm std}_{n,{\rm wht}}^{2} F_{{\rm TIE},0} (A,b) \end{equation} \tag{ 3.1 }$$

where A is a dimensionless parameter depending on the object's optical properties and on the imaging conditions,

$$\begin{equation} A\equiv \pi \gamma \lambda zU_{m}^{2} =(\pi /4)\,\gamma \lambda z/h^{2}, \end{equation} \tag{ 3.2 }$$

and

$$\begin{eqnarray} &&\fl F_{{\rm TIE},0} (A,b)=\int_0^1 {{\rm d}U\int_0^1 {{\rm d}V} \frac{{\rm MTF}_{{\rm det}}^{{\rm 2}} (UU_{m} ,VU_{m} )}{[1+A(U^{2}+V^{2})]^{2}}} . \end{eqnarray} \tag{ 3.3 }$$

Here b is a dimensionless parameter (or a vector of parameters) that describes the MTF of the detector.

It is worth mentioning some useful properties of the function F_TIE,0:

• (1)  
  In the case of an ideal detector, with MTF_det ≡ 1 (b = 0), this function depends only on the dimensionless parameter A and can be very accurately approximated (for any A) by the following formula:

  $$\begin{equation} F_{{\rm TIE},0} (A,0)\approx [1+(4/\pi )A]^{-1}. \end{equation} \tag{ 3.4 }$$

  The relative error of the approximate solution in equation (3.4) does not exceed 1.1%.
• (2)  
  For Gaussian MTF_det, MTF_det(u, v) = exp[−2π² $(\sigma_{\det} /2h)^{2}(u^{2}+v^{2})/U_{m}^{2} ]$, and a contact image, A = 0, the function F_TIE,0 depends only on the dimensionless parameter b = σ_det/(2h) and is written as follows:

  $$\begin{equation} F_{{\rm TIE},0} (0,b)=(\pi /4)[{\rm erf}(2\pi b)/(2\pi b)]^{2}, \end{equation} \tag{ 3.5 }$$

  where erf(x) is the error function, ${\rm erf}(x)\equiv 2/\pi^{1/2}\int_0^x {\exp (-t^{2})\,{\rm d}t}$. For instance, for b = 1/2 equation (3.5) gives F_TIE,0(0, 1/2) ≈ 0.08, and, as expected, F_TIE,0 (0, 0) = 1.
• (3)  
  For large values of A, A ≫ 1, the function F_TIE,0 in equation (3.3) is well approximated by the following simple expression (in which we neglect the effect of the MTF of the detector, compared to the effect of phase retrieval)

  $$\begin{equation} F_{{\rm TIE},0} (A,b)\mathop{{\cong}}_{A\gg 1} F_{{\rm TIE},0} (A,0)\approx (\pi /4)\,A^{-1}. \end{equation} \tag{ 3.6 }$$

In order to quantify the reduction of noise in projection images, as a result of performing phase retrieval, we introduce a gain factor,

$$\begin{eqnarray} &&\fl G_{{\rm TIE},0} (A,b)\equiv \frac{{\rm std}_{n,{\rm cont}}}{{\rm std}_{n,0}}=[F_{{\rm TIE},0} (0,b)/F_{{\rm TIE},0} (A,b)]^{1/2}.\nonumber\\ \end{eqnarray} \tag{ 3.7 }$$

where std_n,0 is the standard deviation of noise in a phase-retrieved phase-contrast image and std_n,cont is the standard deviation of noise in the corresponding unprocessed contact image (of the same object, for the same radiation dose).

Dependencies of the function F_TIE,0 and the corresponding gain factor G_TIE,0 on the parameter A, for various values of parameter b, are shown in figure 1 with solid lines. The dependencies have linear (in the log–log scale) asymptotic at large values of A, namely, F_TIE,0 ∼ A⁻¹ and G_TIE,0 ∼ A^1/2, in accordance with the theoretical results given by equations (3.4)–(3.7). Moreover, in the case of an ideal detector (with MTF_det = 1, which corresponds to b = 0, solid black line in figure 1), equation (3.4) very accurately describes F_TIE,0 and G_TIE,0 for all values of A.

Zoom In Zoom Out Reset image size

TIE is usually valid for λz/h² ⩽ 1 and γ is usually much larger than one, so that typically A is much larger than one, A≫1. For example, for soft materials (composed of light chemical elements) and hard x-rays (with the energies higher than 20 keV), the ratio γ is of the order of 10³. The typical value of σ_det/(2h) for most 2D x-ray detectors is about 1/2. One can see from figure 1 that, in the case of A = 10³, the gain factor is then of the order of 10.

Equation (3.6), which is valid for large values of A, A ≫ 1, allows one to investigate/predict the behaviour of the variance ${\rm std}_{n,0}^{2}$ of noise in a phase-retrieved projection image. Indeed, substituting equation (3.6) in equation (3.1), and taking into account that ${\rm std}_{n,{\rm wht}}^{2} =1/N_{{\rm ph,pxl}}$, one obtains

$$\begin{equation} {\rm std}_{n,0}^{2} \cong \frac{h^{2}}{N_{{\rm ph,pxl}} \gamma \lambda z}. \end{equation} \tag{ 3.8 }$$

• (1)  
  Effect of pixel size, h. If the propagation distance z and x-ray wavelength λ are the independent (primary) parameters describing the imaging system then, according to equation (3.8), the variance of noise is inversely proportional to the photon fluence N_ph,pxl/h² (and, hence, to the radiation dose) and does not depend on the pixel size h. On the other hand, if the Fresnel number N_F = h²/(λz) is the primary parameter describing the imaging system and one fixes it, so that h²/(λz) = const, then, according to equation (3.8), the variance of noise is inversely proportional to the photon fluence and as well as to h². This behaviour of noise is similar to that in the case of conventional radiography. In this case, in order to preserve the level of noise in a contact radiograph (while varying the pixel size of the detector), the incident photon fluence (and radiation dose) should be inversely proportional to h².
• (2)  
  Effect of object-to-detector distance, z. Equation (3.8) indicates that the standard deviation of noise in a phase-retrieved projection image is inversely proportional to the square root of the distance. For instance, in order to decrease the magnitude of noise two times, the distance should be increased four times.
• (3)  
  Energy dependence. There are three terms in equation (3.8) that depend on the x-ray energy, N_ph,pxl, γ and λ. Except for the trivial behaviour of the wavelength, λ ∼ 1/E, both N_ph,pxl and γ depend on object's composition and N_ph,pxl also depends on the thickness of the object. When considering medical imaging, it is the radiation dose delivered to the object that eventually limits the image quality. However, there is always an optimum energy that minimizes the standard deviation of noise in a projection image of an object, for a fixed radiation dose delivered to the object. This optimum energy strongly depends on the imaged object (in terms of its composition as well as its physical dimensions). Note, however, that it is not the noise in an image of an object that is of practical value but rather SNR that quantitatively characterizes the quality of the image (for example, in terms of detectability of some features in the object). This aspect of image quality will be revisited in section 3.5.

Propagation of the variance of random noise from a phase-contrast image to the corresponding contact image obtained with the use of CTF-HOM phase retrieval is described by equation (2.8) with the NPS W_n,0 given by equations (2.18), (2.12) and (2.13),

$$\begin{equation} {\rm std}_{n,0}^{2} \cong {\rm std}_{n,{\rm wht}}^{2} F_{{\rm CTF},0} (a,\gamma ,b), \end{equation} \tag{ 3.9 }$$

where we have introduced a dimensionless parameter a (note that aγ = A, see equation (3.2)),

$$\begin{equation} a\equiv \pi \lambda zU_{m}^{2} =(\pi /4)\,\lambda z/h^{2}, \end{equation} \tag{ 3.10 }$$

and

$$\begin{eqnarray} &&\fl F_{{\rm CTF},0} (a,\gamma ,b)= \int_0^1 {\rm d}U\int_0^1 {\rm d}V \nonumber\\ &&\fl \times\!\frac{\sin^{2}[a(U^{2}+V^{2})+\arctan (\gamma^{-1})]{\rm MTF}_{{\rm det}}^{{\rm 2}} (UU_{m} ,VU_{m} )}{\,(\gamma ^{2}+\!1)\{\sin^{2}[a(U^{2}\!+\!V^{2})+\!\arctan (\gamma^{-1})]+\!\varepsilon (UU_{m},\!VU_{m} )\}^{2}} .\nonumber\\ \end{eqnarray} \tag{ 3.11 }$$

Similar to equation (3.7), a gain factor is introduced as follows:

$$\begin{eqnarray} &&\fl G_{{\rm CTF,}0} (a,\gamma ,b)\equiv \frac{{\rm std}_{n,{\rm cont}}}{{\rm std}_{n,0}}\nonumber\\ &&=[F_{{\rm CTF,0}} (0,\gamma ,b)/F_{{\rm CTF},0} (a,\gamma ,b)]^{1/2}. \end{eqnarray} \tag{ 3.12 }$$

Dependencies of the function F_CTF,0 and the corresponding gain factor G_CTF,0 on the parameter A, for various values of parameter b, are shown in figure 1 with dashed lines. These dependences are conveniently superimposed with the corresponding dependences for the TIE-HOM phase-retrieval algorithm.

One can see from figure 1 that the dependencies for TIE-HOM and CTF-HOM algorithms coincide for relatively small values of A, as long as the condition a ⩽ 1 for the parameter a, defined in equation (3.10), is satisfied. This behaviour is predictable as it is well known [15] that in the near-field regime, which corresponds to large Fresnel numbers, h²/(λz) ⩾ 1, CTF-based reconstruction formula, equation (2.15) converts to the TIE-based reconstruction formula, equation (2.10).

On the other hand, figure 1 indicates that for larger values of a, a > 1, (that is, for example, for longer propagation distances) the function F_CTF,0 as well as the function G_CTF,0, calculated using CTF-HOM algorithm, start to deviate from the corresponding functions F_TIE,0 and G_TIE,0, calculated using TIE-HOM algorithm. Moreover, for each value of parameter γ (which quantifies phase-to-attenuation ratio in the object), there is an asymptotic (maximum) value for the gain factor, (G_CTF,0)_max, and, according to figure 1, this maximum value can be roughly estimated as follows, (G_CTF,0)_max ≈ 0.3γ. For comparison, the TIE-HOM validity condition, h²/(λz) ⩾ 1, starts to be violated at smaller values of the gain factor, about 2–3 times smaller than the corresponding (G_CTF,0)_max.

These results indicate, in particular, that, (1) there might be some benefit of using slightly larger propagation distances than the maximum allowable in TIE, z_TIE,max ≈ h²/λ, and applying CTF-HOM phase-retrieval formula, however, (2) using very large propagation distances, z ≫ z_TIE,max, in combination with CTF-HOM algorithm does not seem to result in any additional benefits, at least in terms of noise reduction/suppression.

Propagation of the variance of random noise from a phase-contrast image to the corresponding contact image obtained with the use of GOA-HOM phase retrieval is described by equation (2.8) with the NPS W_n,0 given by equations (2.26), (2.12) and (2.13),

$$\begin{equation} {\rm std}_{n,0}^{2} \cong {\rm std}_{n,{\rm wht}}^{2} F_{{\rm GOA},0} (c,\omega ,b), \end{equation} \tag{ 3.13 }$$

where

$$\begin{eqnarray} &&\fl F_{{\rm GOA,}0} (c,\omega ,b)=R^{-1}(\omega )\int_0^1 {\rm d}U\int_0^1 {{\rm d}V} \nonumber\\ &&\times\frac{{\rm MTF}_{{\rm det}}^{{\rm 2}} (UU_{m} ,VU_{m} )}{1+c^{2}U^{2}}. \end{eqnarray} \tag{ 3.14 }$$

Here ω is the deviation angle of the analyser crystal from the exact Bragg position, the dimensionless parameter c is defined in equation (2.25) (it plays a role similar to the inverse Fresnel number in PBI in the preceding subsections), and b is, as defined after equation (3.3), a dimensionless parameter (or a vector of parameters) that describes the MTF of the detector. It should be noted that, in order to be consistent with the results of previous subsections, the variance of noise in phase-contrast images, ${\rm std}_{n,{\rm wht}}^{2}$, is supposed to be corrected for the rocking curve value (i.e. $W_{n,{\rm wht}} =h^{2}R(\omega ){\rm std}_{n,{\rm wht}}^{2} )$ and thus it is the variance of noise in the beam before it has been reflected by the analyser crystal.

Similar to equation (3.7), a gain factor is introduced as follows:

$$\begin{eqnarray} &&\fl G_{{\rm GOA},0} (c,\omega ,b)\equiv \frac{{\rm std}_{n,{\rm cont}}}{{\rm std}_{n,0}} \nonumber\\ &&=[R_{0} F_{{\rm GOA},0} (0,0,b)/F_{{\rm GOA},0} (c,\omega ,b)]^{1/2}, \end{eqnarray} \tag{ 3.15 }$$

where std_n,cont is the standard deviation of noise in the contact image (acquired without an analyser crystal). In the case of a Gaussian MTF_det, the function R₀F_GOA,0(0,0,b) is given by equation (3.5).

Regarding the choice of the relative angular spread (Δω/σ) (of diffracted x-rays with respect to the width of the rocking curve of the analyser crystal) in equation (2.25), the following estimations may be indicative. At the x-ray energy of 20 keV (λ = 0.62 × 10⁻⁴ µm), the Si(1 1 1) double-crystal curve has the full-width at half-maximum (FWHM) of 18.341 µrad. The corresponding value for σ is σ = FWHM/1.665 = 11.0 µrad. Assuming h = 1 µm (high resolution detector!) one obtains Δω = λ/(4πh) = 4.93 µrad and hence (Δω/σ) ≈ 0.45. For a bigger pixel size, say h = 10 µm, the same ratio becomes smaller, (Δω/σ) ≈ 0.045.

It should also be noted that at sufficiently high x-ray energies, E > 20 keV, the rocking curve FWHM (and, hence, σ) is proportional to the x-ray wavelength and, therefore, the ratio Δω/σ is energy independent.

Dependencies of the function F_GOA,0 and the corresponding gain factor G_GOA,0 on the parameter ω/σ (relative deviation of the analyser crystal from the Bragg peak position), for various values of the parameter b, γ and Δω/σ, are shown in figure 2.

Zoom In Zoom Out Reset image size

For large values of the parameter c, defined in equation (2.25), the function F_GOA,0(c, ω, b) is well approximated by the following expression:

$$\begin{eqnarray} &&\fl F_{{\rm GOA},0} (c,\omega ,b)\mathop{\cong}_{c\gg 1} R^{-1}(\omega )\int_0^1 {{\rm d}U} \frac{1}{1+c^{2}U^{2}} \nonumber\\ &&\times\int_0^1 {{\rm d}V\,{\rm MTF}_{{\rm det}}^{{\rm 2}} {\rm (}0,VU_{m} )} \approx (\pi /2)R^{-1}(\omega )c^{-1} \nonumber\\ &&\times\int_0^1 {{\rm d}V{\rm MTF}_{{\rm det}}^{{\rm 2}} (0,VU_{m} )} . \end{eqnarray} \tag{ 3.16 }$$

Then, considering the function F_GOA,0 with respect to the argument (ω/σ), it is well approximated by the following expression, {(ω/σ)exp[−(ω/σ)²]}⁻¹, which has its minimum at ω/σ = 2^−1/2 ≈ 0.707. Moreover, for a Gaussian MTF_det, ${\rm MTF}_{{\rm det}} (u,v)=\exp [-2\pi^{2}b^{2}\,(u^{2}+v^{2})/U_{m}^{2} ]$, that depends only on a dimensionless parameter b ≡ σ_det/(2h), the integral in equation (3.16) can be evaluated as follows (see equation (3.5) for details):

$$\begin{eqnarray} &&\fl F_{{\rm GOA},0} (c,\omega ,b)\mathop{\cong}_{c\gg 1} R^{-1}(\omega )c^{-1}(\pi^{3/2}/4){\rm erf}(2\pi b)/(2\pi b).\nonumber\\ \end{eqnarray} \tag{ 3.17 }$$

Considering this last expression for F_GOA,0 as a function of the detector pixel size, h, and assuming that the parameter b is fixed, one can see that F_GOA,0 ∼ h. If, in addition, the photon fluence, N_ph,pxl/h², in a phase-contrast image is fixed (so that the radiation dose delivered to the object is fixed) then, according to equation (2.13), std_n,wht ∼ 1/h and, according to equation (3.13), std_n,0 ∼ h^−1/2. On the other hand, in order to preserve the level of noise in a phase-retrieved projection image, the incident photon fluence (and radiation dose) should be inversely proportional to h.

Using equations (3.15) and (3.17), the gain factor, in the case of large c, can be presented in the following form:

$$\begin{equation} G_{{\rm GOA,0}} (c,\omega ,b)\mathop{\cong}_{c\gg 1} C_{b} [R(\omega )\,c]^{1/2}, \end{equation} \tag{ 3.18 }$$

where the proportionality constant C_b depends only on the detector's parameter(s) b, and in the case of a Gaussian MTF_det (with b = σ_det/(2h)) it is written as follows:

$$\begin{equation} C_{b} =[\pi^{-1/2}{\rm erf}(2\pi b)/(2\pi b)]^{1/2}. \end{equation} \tag{ 3.19 }$$

At b = 0 the constant is C₀ = 0.798, while at b = 1/2 it is C_1/2 = 0.424.

In order to validate the theoretical formalism derived in section 3.1 we performed a series of numerical experiments using a simple numerical phantom. The phantom represents a cylinder with diameter D consisting of a uniform mixture of glandular and adipose tissues [28] in the weight proportion 50/50. X-ray optical properties of the phantom (including the ratio γ = δ/β in table 3) were calculated using x-ray data from the NIST website [29]. Two diameter values have been used in the simulations, D = 2 cm and 8 cm, respectively. In-line phase-contrast images of the cylinder have been simulated by numerically calculating the Fresnel integral applied to the complex amplitude of the wave in the object plane. The propagation distance z was 100 cm. The images were calculated for several x-ray energies including 20, 25, 30, 35, 40 and 45 keV. The projection approximation has been used to calculate the phase shift and the attenuation induced in the transmitted wave by the phantom. The resultant pixel size h in the projections was 10 µm.

A single sinogram has been created for each projection data set. Poisson noise has been generated in the sinograms using certain values for the mean absorbed dose delivered to the phantom. Two doses were used in the numerical simulations, 5 and 500 mGy. For each sinogram, a processed copy was created using a 1D implementation of TIE-HOM phase-retrieval algorithm (so that a single row of the detector was phase retrieved for each projection view). The resultant phase-retrieved projections (2k and 8k pixel wide for 2 cm and 8 cm phantoms, respectively) contained intensity values and have been analysed by calculating the standard deviation in the centre of the projections (where the attenuation was the strongest) and in the background (where the attenuation was zero). The results of this analysis are summarized in tables 1 and 2.

Table 1. Standard deviation of noise in phase-contrast (std_n,wht) and phase-retrieved (std_n,0) projection images of the cylindrical phantom with the diameter of 2 cm. Superscripts (e) and (c) correspond to pixels outside (exterior) and in the centre of the cylinder's projection, respectively.

	E (keV)	20	25	30	35	40	45
Dose	I(e)	1	1	1	1	1	1
	I(c)	0.2742	0.4281	0.5215	0.5777	0.6131	0.6368
500 mGy	${\rm std}_{n,{\rm wht}}^{({\rm e})}$	0.037 32	0.032 02	0.027 78	0.024 65	0.022 59	0.021 26
	${\rm std}_{n,{\rm wht}}^{({\rm c})}$	0.019 54	0.020 95	0.020 06	0.018 74	0.017 69	0.016 97
	${\rm std}_{n,0}^{({\rm e})}$	0.005 966	0.005 571	0.004 896	0.004 331	0.004 083	0.003 998
	${\rm std}_{n,0}^{({\rm c})}$	0.003 195	0.003 465	0.003 394	0.003 266	0.003 241	0.003 193
5 mGy	${\rm std}_{n,{\rm wht}}^{({\rm e})}$	0.3732	0.3202	0.2778	0.2465	0.2259	0.2126
	${\rm std}_{n,{\rm wht}}^{({\rm c})}$	0.1954	0.2095	0.2006	0.1874	0.1769	0.1697
	${\rm std}_{n,0}^{({\rm e})}$	0.062 49	0.053 03	0.047 04	0.044 13	0.041 72	0.039 20
	${\rm std}_{n,0}^{({\rm c})}$	0.031 72	0.034 02	0.034 04	0.033 23	0.032 39	0.032 09

Table 2. Standard deviation of noise in phase-contrast (std_n,wht) and phase-retrieved (std_n,0) projection images of the cylindrical phantom with the diameter of 8 cm. Superscripts (e) and (c) correspond to pixels outside (exterior) and in the centre of the cylinder's projection, respectively.

	E (keV)	20	25	30	35	40	45
Dose	I(e)	1	1	1	1	1	1
	I(c)	0.005 653	0.033 61	0.073 98	0.1114	0.1413	0.1645
500 mGy	${\rm std}_{n,{\rm wht}}^{({\rm e})}$	0.050 56	0.051 21	0.049 51	0.047 13	0.044 98	0.043 26
	${\rm std}_{n,{\rm wht}}^{({\rm c})}$	0.003 801	0.009 388	0.013 47	0.015 73	0.016 91	0.017 54
	${\rm std}_{n,0}^{({\rm e})}$	0.008 114	0.008 324	0.008 404	0.008 298	0.008 133	0.008 141
	${\rm std}_{n,0}^{({\rm c})}$	0.000 6213	0.001 548	0.002 279	0.002 752	0.003 075	0.003 312
5 mGy	${\rm std}_{n,{\rm wht}}^{({\rm e})}$	0.5056	0.5121	0.4951	0.4713	0.4498	0.4326
	${\rm std}_{n,{\rm wht}}^{({\rm c})}$	0.038 01	0.093 88	0.1347	0.1573	0.1691	0.1754
	${\rm std}_{n,0}^{({\rm e})}$	0.082 08	0.084 41	0.082 79	0.082 16	0.082 57	0.081 55
	${\rm std}_{n,0}^{({\rm c})}$	0.006 249	0.015 44	0.022 72	0.027 54	0.030 77	0.033 05

Figure 3 shows experimental gain factors calculated using standard deviations from tables 1 and 2 as well as theoretical values (solid black line) calculated assuming 1D phase retrieval with TIE-HOM. It should be emphasized that the gain factors calculated using the results of the numerical experiments are in good agreement with those obtained using our theoretical formulae from section 3.1. Moreover, though the standard deviation of noise varies across a phase-contrast image and the corresponding phase-retrieved image (due to the average intensity variation in the images), the gain factor is shift invariant. This is qualitatively explained below.

Zoom In Zoom Out Reset image size

As was mentioned in section 2.1, if the intensity in an image is slowly varying (as in most real images) then the corresponding noise in the image is quasi-stationary [9], with the autocorrelation function well approximated by a product of two factors: a slowly varying factor due to slow variations in average intensity and a short-range function describing correlation between neighbouring points. Phase retrieval using TIE-HOM preserves this quasi-stationarity because its point-spread function, P_TIE(x, y) (corresponding to the Fourier filter in equation (2.10)), is essentially spatially localized:

$$\begin{equation*} P_{{\rm TIE}} (x,y)=P_{{\rm TIE}} (r)=K_{0} (r/l_{{\rm TIE}} )/(2\pi l_{{\rm TIE}}^{2} ), \end{equation*}$$

where K₀ is the zeroth order modified Bessel function of the second kind, r = (x² + y²)^1/2 is the radial coordinate and l_TIE = (1/2)(γλz/π)^1/2 is the characteristic length scale of the variation of this point-spread function (note that function K₀ exponentially decreases at large values of its argument and its standard deviation is equal to 2l_TIE).

Although we present numerical results only for propagation-based phase-contrast images reconstructed with TIE-HOM algorithm, similar quasi-stationarity of noise is expected in the case of analyser-based phase-contrast images reconstructed with the GOA-HOM algorithm, equation (2.23). The point-spread function P_GOA(x) of the GOA-HOM imaging method, corresponding to the Fourier filter in equation (2.23), is written as follows:

$$\begin{equation*} P_{{\rm GOA}} (x)=\left\{{\begin{array}{@{}ll@{}} l_{{\rm GOA}}^{-1} \exp (-x/l_{{\rm GOA}} ), & x>0, \\ 0, & x<0, \\ \end{array}} \right. \end{equation*}$$

where l_GOA = γ(ω/σ) (kσ)⁻¹ is the characteristic length scale of variation of this point-spread function. Note that (kσ)⁻¹ is essentially the projection of the extinction depth of the analyser crystal onto the object plane.

It should be emphasized that knowledge of standard deviation of noise in an image is not sufficient in order to quantify image quality. A relative characteristic, such as, for example, SNR is required for this purpose. In our analysis below, we use two definitions for SNR in projection images. In the first definition, mean intensity in a pixel is treated as signal and the standard deviation of noise in the same pixel is treated as noise, so that the SNR is simply, SNR = I/std_n.

Alternatively, SNR can be defined as follows:

$$\begin{equation*} {\rm SNR'}=I\mu t/{\rm std}_{n} , \end{equation*}$$

where Iμt is the contrast in the projection image associated with a small feature (of thickness t) located inside or on the surface of the object and μ is either the difference of the linear attenuation coefficient of the feature and of the object or the linear attenuation coefficient of the feature, respectively.

Using equation (3.8), one can write, in the case of TIE-HOM (and large A, see table 3)

$$\begin{equation} {\rm SNR}_{{\rm TIE}} \cong I^{3/2}h^{-1}(\gamma \lambda zN_{{\rm inc,pxl}} )^{1/2}. \end{equation} \tag{ 3.20 }$$

Similarly, using equations (3.13) and (3.17) and assuming perfect detector (so that b = 0 in equation (3.17)), one has

$$\begin{equation} {\rm SNR}_{{\rm GOA}} \cong I^{3/2}[(2/\pi )R(\omega )cN_{{\rm inc,pxl}} ]^{1/2}, \end{equation} \tag{ 3.21 }$$

where, according to equation (2.25), c = 2π γ(ω/σ) (Δω/σ).

Using the optimum deviation angle, ω_opt = 2^−1/2σ, that minimizes the standard deviation of noise and hence maximizes SNR, and assuming Si(111) reflection in the monochromator as well as in the analyser crystal, the product R(ω)c in equation (3.21) is written as follows, 1.2076γ/h[µm].

It should be noted that in the case of phase retrieval using TIE-HOM (with A ≫ 1) SNR_TIE does not depend on the pixel size h (as soon as propagation distance z is fixed), while in the case of phase retrieval using GOA-HOM (with b = const) SNR_GOA is proportional to h^1/2. Note, however, that in the case of phase retrieval using TIE-HOM and fixing the Fresnel number, N_F = h²/(λ z), instead of propagation distance z, SNR_TIE is proportional to the pixel size, exactly as in conventional (contact) images (see table 4 for details).

Table 4. Comparison of several key characteristics of noise in a phase-retrieved image, for three imaging methods (including conventional, PBI and ABI). Here Φ_ph,inc is the photon fluence (the number of photons incident per unit cross-sectional area) in the incident beam. Realistic detector with b ≡ σ_det/(2h) = 1/2 is assumed.

Method	Conventional	PBI: TIE-HOMa	ABI: GOA-HOMb
${\rm std}_{n,0}^{2}$	∼1/(Φph,inc × h2)	∼1/Φph,inc	∼1/(Φph,inc × h)
SNR2	∼(Φph,inc × h2)	∼Φph,inc	∼(Φph,inc × h)
Gain for γ = 103	1	10	≈16
Gain for γ = 102	1	≈3.5	≈5

^aAssuming fixed z. If the Fresnel number is fixed then the h-dependence coincides with that of the conventional CT reconstruction. Gain values, $G_{{\rm TIE,0}} (A\gg 1,b=1/2)\cong \sqrt {0.08\gamma \lambda z/h^{2}}$, correspond to a ≡ (π/4)λz/h² = 1 (so that A = γ). ^bGain values, G_GOA,0(c ≫ 1, b = 1/2) ≅ 0.424(Rc)^1/2, where c = 2πγ(ω/σ)(Δω/σ), correspond to Δω/σ = 1/2 (high resolution detector, h ≈ 1 µm) and ω/σ = (ω/σ)_opt = 2^−1/2 (with the corresponding R = exp(−1/2)).

For a given mean absorbed dose D_abs delivered to an object, the corresponding number of photons N_inc,pxl incident onto a normal surface area S = h² (equal to the area of one detector pixel) of the object is calculated as follows [30]:

$$\begin{eqnarray} &&\fl N_{{\rm inc,pxl}} [{\rm ph}\,{\rm pxl}^{-1}]=D_{{\rm abs}} [{\rm mGy}]\frac{h^{2}[{\rm cm}]L[{\rm cm}]\rho [{\rm g}\,{\rm cm}^{-3}]}{1.602\times 10^{-10}E[{\rm keV}]R_{{\rm abs}}},\nonumber\\ \end{eqnarray} \tag{ 3.22 }$$

where R_abs is the fraction of the incident energy absorbed by the object, L and ρ is the thickness and the mass density of the object, respectively.

We applied equations (3.20)–(3.22) to the same numerical phantom that has been used in section 3.4. Figure 4 shows SNR versus x-ray energy, for two values of the phantom's diameter, D = 2 and 8 cm, for the radiation dose D_abs = 1 mGy. Figures 4(a), (c), (e) and (g) show SNR values in the centre of projection images (where the attenuation is the strongest) calculated using equations (3.20) and (3.21), while figures 4(b), (d), (f) and (h) show SNR values calculated using the alternative definition, SNR' = SNR × (μt), assuming that a small feature (of thickness t = 0.5 cm) is located on the surface of the cylinder and has the same linear attenuation coefficient μ as the cylinder.

Zoom In Zoom Out Reset image size

Apparently, there is an optimum energy (which depends on the object's attenuation properties and size as well as on the imaging method and the corresponding phase-retrieval algorithm) at which SNR is maximal. Figure 4 indicates that (i) the optimum energy increases with the object's size; (ii) the corresponding maximal SNR quickly decreases with the object's size; (iii) the optimum energy that maximizes SNR is slightly higher compared to the optimum energy that maximizes SNR'; (iv) the optimum energy that maximizes SNR in an image reconstructed using GOA-HOM algorithm is either higher compared to that for TIE-HOM phase retrieval (if the propagation distance is fixed, figures 4(c) and (d)) or coincides with that for TIE-HOM phase retrieval (if the Fresnel number is fixed, figures 4(g) and (h)).

Figure 4 also shows that, for the chosen detector pixel size, h = 10 µm, and free-space propagation distance (used in TIE-HOM), z = 100 cm, propagation-based image reconstructed using TIE-HOM produces slightly larger maximal SNR than the corresponding analyser-based image reconstructed using GOA-HOM, especially for the smaller object. This difference becomes even larger if, instead of fixing the propagation distance (as in figures 4(c) and (d)), one fixes the Fresnel number, N_F = 1 (as in figures 4(g) and (h)).

Phase retrieval improves (increases) SNR (in comparison with the corresponding contact image) by the same gain factor G₀ that has been introduced above, in sections 3.1–3.3, in order to quantify the improvement (decrease) of the standard deviation of noise. According to table 3 (see the rows corresponding to the TIE-HOM phase retrieval, with 100 cm propagation distance and 10 µm pixel size) the gain factor is of the order of 5–6 in the case of 1D phase retrieval and 25–33 in the case of 2D phase retrieval (shown in figure 4). For comparison, in the case of the GOA-HOM phase retrieval (with 10 µm pixel size), the corresponding gain factor is in the range 12–13.

Importantly, unlike the standard deviation of noise, the gain factor is dose independent. It depends only on the phase-retrieval algorithm (through its Fourier filter) and detector's properties (via detector's MTF).

Also, we would like to emphasize the fundamental difference, in terms of the relationship between the variance of noise and the photon statistics, of the three imaging methods, including conventional radiography, PBI combined with TIE-HOM phase retrieval and ABI combined with GOA-HOM phase retrieval. These relationships follow from our theoretical analysis presented in sections 3.1 and 3.3 as well as in this section, and are summarized in table 4. Also shown in table 4 are the typical gain values for two PCI methods.

Using equations (2.31)–(2.32) and (2.7) and ramp filter, G(ξ) = ξ, one easily obtains the following analytical expression for the variance of noise in a CT-reconstructed slice containing the spatial distribution of the linear attenuation coefficient, assuming uncorrelated (white) Poisson noise in contact projection images with the standard deviation std_n,wht:

$$\begin{equation} {\rm std}_{n,\mu ,{\rm NN}}^{2} =\frac{\pi^{2}/12}{N_{{\rm a}} h^{2}}{\rm std}_{n,{\rm wht}}^{2} \end{equation} \tag{ 4.1 }$$

Equation (4.1) was found to be in excellent agreement with the numerical results obtained using an implementation of the FBP CT reconstruction algorithm with NN interpolation (in the backprojection step) and applied to flat-field (i.e. without an object) projections with generated Poisson noise of different levels (in the range of 0.5–4%).

More specifically, we used 900 projections 1k×1k pixels each, with the pixel size h = 1 µm, and reconstructed a slice containing distribution of the imaginary part β of the complex refractive index at the wavelength λ = 10⁻⁴ µm [μ = (4π/λ)β]. Using the central part of the reconstructed slice, 256 × 256 pixels, we calculated the standard deviation of the reconstruction noise, for five levels of Poisson noise in the input sinogram/projections. The results of this numerical analysis are plotted as circles in figure 5 together with the result of linear fit (solid line). From this fit, the proportionality constant in equation (4.1) has been estimated² to be C = 0.816 ± 0.005 which perfectly coincides with the theoretical value, π²/12 = 0.822.

Zoom In Zoom Out Reset image size

It is worth mentioning that an implementation of the FBP algorithm with LI in the backprojection step produces less noise in the reconstructed slices so that the proportionality constant C in similar numerical simulations was significantly smaller, namely, C = 0.382 ± 0.002. This result is in agreement with a theoretical formula for the variance of noise in a reconstructed slice, obtained in [27] using discrete representation of the ramp filtering:

$$\begin{equation} {\rm std}_{n,\mu ,{\rm LI}}^{2} =\frac{\pi^{2}-3}{18N_{{\rm a}} h^{2}}{\rm std}_{n{\rm ,wht}}^{2} . \end{equation} \tag{ 4.2 }$$

Equation (4.2) can be easily obtained using equations (2.31), (2.33) and (2.7) with ramp filter, G(ξ) = ξ.

For our further analysis we introduce a NPS W_n,μ(u, v, w) for 3D noise distribution in the CT-reconstructed volume containing the linear attenuation coefficient. This 3D NPS is related to 2D NPS W_n,0(ξ, w) in the corresponding contact 2D projections as follows (hereafter we restrict our consideration to the case of ramp filter and NN interpolation),

$$\begin{eqnarray} &&\fl W_{n,\mu} (u,v,w)=\frac{\pi}{N_{{\rm a}}}\xi W_{n,0} (\xi ,w),\tqs \xi =\sqrt {u^{2}+v^{2}} .\nonumber\\ \end{eqnarray} \tag{ 4.3 }$$

Then the variance of noise in the reconstructed volume is written as follows (compare with equation (2.31)):

$$\begin{equation} {\rm std}_{n,\mu}^{2} =\frac{2\pi^{2}}{N_{{\rm a}}}\int_0^{U_{m}} {{\rm d}\xi \xi^{2}\int_{-U_{m}}^{+U_{m}} {{\rm d}wW_{n,0} (\xi ,w)}} . \end{equation} \tag{ 4.4 }$$

Below we investigate the effect of the three phase-retrieval algorithms described in section 2 on noise propagation from phase-contrast projections to the reconstructed volume.

Using equations (4.4), (2.11)–(2.13), the variance of noise in a reconstructed volume is related to the variance of uncorrelated noise in the corresponding phase-contrast projection images as follows:

$$\begin{equation} {\rm std}_{n,\mu}^{2} ={\rm std}_{n,{\rm wht}}^{2} \frac{\pi ^{2}/12}{N_{{\rm a}} h^{2}}F_{{\rm TIE}} (A,b), \end{equation} \tag{ 4.5 }$$

where A and b are defined in section 3.1 (see equation (3.2) and succeeding text) and

$$\begin{eqnarray} &&\fl F_{{\rm TIE}} (A,b)=3\int_0^1 {{\rm d}\Xi \Xi^{2}\int_0^1 {{\rm d}W\frac{{\rm MTF}_{{\rm det}}^{{\rm 2}} (\Xi U_{m} ,WU_{m} )}{[1+A(\Xi ^{2}+W^{2})]^{2}}}} .\nonumber\\ \end{eqnarray} \tag{ 4.6 }$$

Equation (4.5) is a direct extension of equation (4.1) to the case of CT reconstruction from phase-contrast projection data; it takes into account the resolution of the detector (via its MTF) as well as the TIE-HOM phase retrieval.

Below we provide some analysis on the function F_TIE, similar to that carried out in section 3.1:

• (1)  
  In the case of a perfect detector, with MTF_det ≡ 1 (b = 0), the function F_TIE depends only on the dimensionless parameter A. Moreover, for large values of A, i.e. A ≫ 1, F_TIE can be approximated as follows:

  $$\begin{equation} F_{{\rm TIE}} (A\gg 1,0)\cong (3\pi /8)\frac{\ln A-1}{A^{2}}. \end{equation} \tag{ 4.7 }$$
• (2)  
  For a Gaussian MTF_det, ${\rm MTF}_{{\rm det}} (u,v)=\exp [-2\pi^{2}b^{2}(u^{2}+v^{2})/U_{m}^{2} ]$, and a contact image, A = 0, the function F_TIE depends only on the dimensionless parameter b ≡ σ_det/(2h) and is written as follows:

  $$\begin{eqnarray} &&\fl F_{{\rm TIE}} (0,b)=\frac{3\sqrt \pi}{2}\frac{{\rm erf}(2\pi b)}{2\pi b} \nonumber\\ &&\times\left\{{\frac{\sqrt \pi}{4}\frac{{\rm erf}(2\pi b)}{(2\pi b)^{3}}-\frac{\exp [-(2\pi b)^{2}]}{2(2\pi b)^{2}}} \right\}. \end{eqnarray} \tag{ 4.8 }$$

  For instance, for b = 1/2 equation (4.8) gives F_TIE (0, 1/2) ≈ 0.012 (for comparison, F_TIE (0, 0) = 1).
• (3)  
  For large values of A, i.e. A ≫ 1, and Gaussian MTF_det, the function F_TIE in equation (4.6) is well approximated by the following expression (that unlike equation (3.6) depends not only on A but also on b):

  $$\begin{eqnarray} &&\fl F_{{\rm TIE}} (A,b)\mathop{\cong}_{A\gg 1} (3\pi /8) \nonumber\\ &&\times \frac{\ln (A)-1-\{\gamma _{E} +\ln [(2\pi b)^{2}]+Ei[1,(2\pi b)^{2}]\}}{A^{2}},\nonumber\\ \end{eqnarray} \tag{ 4.9 }$$

where γ_E = 0.577 216 is the Euler's constant and Ei(n, x) is the exponential integral, $Ei(n,x)=\int_1^\infty {{\rm d}t\,{\rm e}^{-xt}t^{-n}}$. Considering the sum in the curly brackets, γ_E + ln(x) + Ei(1, x), it is worth mentioning that (i) for large values of x, x > 2 (b ⩾ 1/4), the exponential integral is small compared to the first two terms and can be neglected, so that γ_E + ln(x) + Ei(1, x) ≈ γ_E + ln(x); and (ii) in the opposite case of small x, x < 0.1 (b ⩽ 1/20), this sum is well approximated as follows: γ_E + ln(x) + Ei(1, x) ≈ x. Using this latter property, one can see that in the case of ideal detector (with b = 0), equation (4.9) converts to equation (4.7).

In order to quantify the reduction of noise in CT-reconstructed slices as a result of applying TIE-HOM phase retrieval to phase-contrast projections, prior to FBP CT reconstruction (with respect to CT reconstruction from unprocessed contact projection images of the same object with the same radiation dose), we introduce the following gain factor:

$$\begin{eqnarray} &&\fl G_{{\rm TIE}} (A,b)\equiv \frac{{\rm std}_{n,\mu} (0,b)}{{\rm std}_{n,\mu} (A,b)}=[F_{{\rm TIE}} (0,b)/F_{{\rm TIE}} (A,b)]^{1/2}.\nonumber\\ \end{eqnarray} \tag{ 4.10 }$$

Dependencies of the function F_TIE and the corresponding gain factor G_TIE on the parameter A, for various values of the parameter b (describing Gaussian MTF of the detector), calculated using equations (4.6) and (4.10), are shown in figure 6 with solid lines. The dependencies have almost linear (in the log–log scale) asymptotics at large values of A, in accordance with the theoretical results given in equations (4.8)–(4.10).

Zoom In Zoom Out Reset image size

TIE is usually valid for λz/h² ⩽ 1 and γ is usually much bigger than one, so that A ≫ 1.

For example, for soft materials (composed of light elements) and hard x-rays (with the energy larger than about 20 keV), the value of γ is typically of the order of 10³. The typical value of the ratio σ_det/(2h) is about 1/2. One can see from figure 6 that, in the case of A = 10³, the gain factor is of the order of 50–60.

Equation (4.9), valid for large values of A(A ≫ 1), allows one to present the variance of noise in a reconstructed slice as follows:

$$\begin{eqnarray} &&\fl {\rm std}_{n,\mu}^{2} \cong (\pi /2)\frac{\ln [(\pi /4)\gamma \lambda z/h^{2}]-1-C(b)}{N_{{\rm a}} (\gamma \lambda z)^{2}}\frac{h^{2}}{N_{{\rm ph,pxl}}},\nonumber\\ \end{eqnarray} \tag{ 4.11 }$$

where C(b) = γ_E + ln[(2πb)²] + Ei[1, (2πb)²] is a constant that depends only on the detector's parameter b.

• (1)  
  Effect of the pixel size. If the photon fluence, N_ph,pxl/h², in the phase-contrast images is fixed (so that the radiation dose delivered to the object is fixed) then, according to equation (4.11), the standard deviation of noise in the reconstructed slice depends on the pixel size via the following term: std_n,μ ∼ [ln(A) − 1 − C(b)]^1/2. For instance, assuming b = 1/2 and hence C = 2.8667, then if the original value of A was 10³ and the pixel size is increased by a factor of 2 so that new value of A becomes 250 then the standard deviation std_n,μ would be reduced by a factor of [(ln 10³ − 3.8667)/(ln 250 − 3.8667)]^1/2 ≈ 1.3556. For comparison, in the case of a perfect detector, with b = 0 and hence C = 0, the corresponding factor would be [(ln 10³ − 1)/(ln 250 − 1)]^1/2 ≈ 1.1431. Thus, in a small range of large A (≫1), the standard deviation of noise is almost independent of the pixel size.
• (2)  
  Effect of object-to-detector distance, z. According to equation (4.11), the standard deviation of noise in a reconstructed slice is approximately inversely proportional to the distance.
• (3)  
  Energy dependence. There are three variables in equation (4.11) that depend on the x-ray energy, N_ph,pxl, γ and λ. The same qualitative analysis, which has been presented in section 3.1, is applicable to this case as well. More detailed analysis of the energy dependence of SNR will be presented in section 4.6.

Using equations (4.4), (2.18), (2.12)–(2.13), the variance of noise in a reconstructed volume is related to the variance of Poisson noise in the corresponding phase-contrast projection images as follows:

$$\begin{equation} {\rm std}_{n,\mu}^{2} \cong {\rm std}_{n,{\rm wht}}^{2} \frac{\pi ^{2}/12}{N_{{\rm a}} h^{2}}F_{{\rm CTF}} (a,\gamma ,b), \end{equation} \tag{ 4.12 }$$

where a is defined in equation (3.10) and, as above, b is a dimensionless parameter that describes the MTF of the detector,

$$\begin{eqnarray} &&\fl F_{{\rm CTF}} (a,\gamma ,b)= 3\int_0^1 {\rm d}\Xi \,\Xi^{2}\int_0^1 {\rm d}W \sin^{2}[a(\Xi ^{2}+W^{2}) \nonumber\\ &&+\arctan (\gamma^{-1})]{\rm MTF}_{{\rm det}}^{{\rm 2}} (\Xi U_{m} ,WU_{m} ) \nonumber\\ &&\times\,{[}(\gamma^{2}+1)\{\sin^{2}[a(\Xi^{2}+W^{2})+\arctan (\gamma^{-1})] \nonumber\\ &&+\,\varepsilon (\Xi U_{m} ,WU_{m} )\}^{2}]^{-1} , \end{eqnarray} \tag{ 4.13 }$$

where the regularization parameter ε is defined in equation (2.17).

Similar to equation (4.10), we introduce a gain factor,

$$\begin{eqnarray} &&\fl G_{{\rm CTF}} (a,\gamma,b)\equiv \frac{{\rm std}_{n,\mu} (0,\gamma ,b)}{{\rm std}_{n,\mu} (a,\gamma ,b)} \nonumber\\ &&=[F_{{\rm CTF}} (0,\gamma ,b)/F_{{\rm CTF}} (a,\gamma ,b)]^{1/2}. \end{eqnarray} \tag{ 4.14 }$$

Dependencies of the function F_CTF and the corresponding gain function G_CTF on the parameter A = aγ, for γ = 10² and 10³ and various values of parameter b of the Gaussian MTF of the detector, calculated using equations (4.13) and (4.14), are shown in figure 6 with dashed lines. These dependences are superimposed with the corresponding dependencies for a combination of the TIE-HOM phase-retrieval algorithm with the FBP CT reconstruction.

The behaviour of the function F, as well as the gain factor G, is similar to that observed in the case of phase retrieval using CTF-HOM and investigated in section 3.2. Indeed, the curves in figure 6, corresponding to TIE-HOM and CTF-HOM algorithm, coincide for relatively small values of A, as long as parameter a is small enough, a ⩽ 1. Similarly, for larger values of a, a > 1, the function F_CTF, as well as the gain factor G_CTF, calculated using CTF-HOM algorithm, starts to deviate from the corresponding function F_TIE and G_TIE, respectively, calculated using TIE-HOM algorithm. Moreover, for each value of parameter γ (which quantifies phase-to-attenuation ratio in the object), there is an asymptotic ('maximum') value for the gain factor, (G_CTF)_max, and, according to figure 6, this 'maximum' value can be roughly estimated as follows: (G_CTF)_max ≈ 0.3γ (note also that the TIE validity condition, h²/(λz) > 1, usually starts to be violated at values of the gain factor smaller than 0.3γ). This value represents the largest possible gain factor, i.e. the reduction of noise level, that can be achieved in the reconstructed CT data obtained using PBI with CTF-HOM phase retrieval in combination with FBP CT reconstruction, compared to conventional CT reconstruction from contact projection images.

The variance of noise in a reconstructed μ-slice is related to the variance of white noise in the corresponding phase-contrast images as follows:

$$\begin{equation} {\rm std}_{n,\mu}^{2} \cong {\rm std}_{n,{\rm wht}}^{2} \frac{\pi ^{2}/12}{N_{{\rm a}} h^{2}}F_{{\rm GOA}} (c,\omega ,b), \end{equation} \tag{ 4.15 }$$

where c is defined previously, in equation (2.25).

The corresponding gain factor is defined as follows:

$$\begin{eqnarray} &&\fl G_{{\rm GOA}} (c,\omega ,b)\equiv \frac{{\rm std}_{n,\mu ,{\rm cont}}}{{\rm std}_{n,\mu} (c,\omega ,b)} \nonumber\\ &&=[R_{0} F_{{\rm GOA}} (0,0,b)/F_{{\rm GOA}} (c,\omega ,b)]^{1/2}, \end{eqnarray} \tag{ 4.16 }$$

where std_n,μ,cont corresponds to a contact image (acquired without an analyser crystal). In the case of a Gaussian MTF_det, the function R₀F_GOA(0,0,b) coincides with that given by equation (4.8).

Regarding the function F_GOA, we consider the following two cases which are different in terms of the relative orientation of the CT rotation axis and the diffraction plane of the analyser crystal.

• (1)  
  Rotation axis is parallel to the diffraction plane of the analyser crystal:

  $$\begin{eqnarray} &&\fl F_{{\rm GOA},\vert \vert} (c,\omega ,b)=3R^{-1}(\omega )\int_0^1 {\rm d}\Xi \,\Xi^{2}\int_0^1 {\rm d}W \nonumber\\ &&\times\frac{\,{\rm MTF}_{{\rm det}}^{{\rm 2}} (\Xi U_{m} ,WU_{m} )}{1+c^{2}W^{2}}. \end{eqnarray} \tag{ 4.17 }$$

  Recall that according to equation (2.25) c = 2π γ(ω/σ) (Δω/σ), where Δω ≡ k⁻¹U_m = λ/(4πh). Assuming that c ≫ 1, the inner integral in equation (4.17) can be simplified:

  $$\begin{eqnarray*}&&\fl F_{{\rm GOA,}\parallel} (c,\omega ,b)\mathop{\cong}_{c\gg 1} 3R^{-1}(\omega ) \nonumber\\ &&\times\int_0^1 {{\rm d}\Xi \,\Xi^{2}{\rm MTF}_{{\rm det}}^{{\rm 2}} (\Xi U_{m} ,0)\int_0^1 {{\rm d}W\frac{1}{1+c^{2}W^{2}}}} \end{eqnarray*}$$

  and it can be evaluated as follows:

  $$\begin{eqnarray} &&\fl F_{{\rm GOA},\vert \vert} (c,\omega ,b)\mathop{\cong}_{c\gg 1} (3\pi /2)R^{-1}(\omega )c^{-1} \nonumber\\ &&\times \int_0^1 {{\rm d}\Xi \,\Xi^{2}{\rm MTF}_{{\rm det}}^{{\rm 2}} (\Xi U_{m} ,0)} . \end{eqnarray} \tag{ 4.18 }$$

  Considering the function F_GOA,∥ in equation (4.18) with respect to its argument (ω/σ), it is easy to show that it has its minimum at the following deviation angle ω_opt,|| of the analyser crystal,

  $$\begin{equation} \omega_{{\rm opt},\vert \vert} =\sigma /\sqrt 2 . \end{equation} \tag{ 4.19 }$$

  This deviation angle corresponds to the rocking curve value R_opt,∥ = R₀ exp(−1/2) = 0.6065R₀. It is worth noting that the corresponding gain factor, G_GOA,||(c, ω, b), has its maximum at the deviation angle ω_opt,|| of the analyser crystal.For a Gaussian MTF of the detector, equation (4.18) can be evaluated analytically as follows:

  $$\begin{eqnarray*}&&\fl F_{{\rm GOA},\vert \vert} (c,\omega ,b)\mathop{\cong}_{c\gg 1} R^{-1}(\omega )c^{-1}(3\pi /2) \nonumber\\ &&\times\left\{{\frac{\sqrt \pi}{4}\frac{{\rm erf}(2\pi b)}{(2\pi b)^{3}}-\frac{\exp [-(2\pi b)^{2}]}{2(2\pi b)^{2}}} \right\}. \end{eqnarray*}$$

  Also, assuming that b is fixed, the function F_GOA,∥ depends on the pixel size h as follows: F_GOA,|| ∼ h. Therefore, for a fixed photon fluence in the image, N_ph,pxl/h², so that std_n,wht ∼ 1/h, the standard deviation of noise in the reconstructed slice is related to the pixel size as follows:

  $$\begin{equation} {\rm std}_{n,\mu ,\vert \vert} =\frac{{\rm const}}{\sqrt {N_{{\rm a}}} h}\frac{{\rm std}_{n,{\rm wht}}}{\sqrt {F_{{\rm GOA},\vert \vert}}}\sim \frac{1}{h^{3/2}}. \end{equation} \tag{ 4.20 }$$

  Equation (4.20) indicates that, despite the fact that the gain factor (due to the phase retrieval as opposed to a pure absorption image) degrades with the pixel size (as shown previously in section 3.3), noise in the reconstructed slice decreases with the pixel size (due to better photon statistics in an individual detector's pixel, the std_n,wht term, and due to CT reconstruction, the first factor in equation (4.20)).Regarding the corresponding gain function at large values of c, the following result is useful for simple quantitative assessments:

  $$\begin{equation} G_{{\rm GOA},\vert \vert} (c,\omega ,b)\mathop{\cong}_{c\gg 1} C_{b,\vert \vert} R^{1/2}(\omega )c^{1/2}, \end{equation} \tag{ 4.21 }$$

  where the proportionality constant C_b,∥ depends only on the detector's parameter b, and in the case of a Gaussian MTF of the detector (with b = σ_det/(2h)) it is written as follows:

  $$\begin{equation*} C_{b,\vert \vert} =[\pi^{-1/2}{\rm erf}(2\pi b)/(2\pi b)]^{1/2}. \end{equation*}$$

  At b = 0 the constant is C_0,∥ = 0.798, while at b = 1/2 it is C_1/2,|| = 0.424.
• (2)  
  Rotation axis is perpendicular to the diffraction plane of the analyser crystal:

  $$\begin{eqnarray} &&\fl F_{{\rm GOA},\bot} (c,\omega ,b)=3R^{-1}(\omega )\int_0^1 {\rm d}\Xi \frac{\Xi^{2}}{1+c^{2}\Xi^{2}} \nonumber\\ &&\times\int_0^1 {{\rm d}W\,{\rm MTF}_{{\rm det}}^{{\rm 2}} (\Xi U_{m} ,WU_{m} )} . \end{eqnarray} \tag{ 4.22 }$$

Assuming again that c ≫ 1, the function F_GOA,⊥ can accurately be approximated as follows³:

$$\begin{eqnarray} &&\fl F_{{\rm GOA},\bot} (c,\omega ,b)\mathop{\cong}_{c\gg 1} 3R^{-1}(\omega )c^{-2}\nonumber\\ &&\times\int_0^1 {{\rm d}\Xi \int_0^1 {{\rm d}W{\rm MTF}_{{\rm det}}^{{\rm 2}} (\Xi U_{m} ,WU_{m} )}} . \end{eqnarray} \tag{ 4.23 }$$

Then, considering the function F_GOA,⊥ with respect to its argument (ω/σ), it is easy to show that it has its minimum at the following deviation angle ω_opt,⊥ of the analyser crystal:

$$\begin{equation} \omega_{{\rm opt},\bot} \approx \sigma . \end{equation} \tag{ 4.24 }$$

This deviation angle corresponds to the rocking curve value R_opt,⊥ = R₀ exp(−1) = 0.3679R₀. It is worth noting that the corresponding gain factor, G_GOA,⊥(c, ω, b), has its maximum at the deviation angle ω_opt,⊥ of the analyser crystal.

For a Gaussian MTF of the detector, equation (4.23) can be evaluated analytically as follows:

$$\begin{eqnarray} &&\fl F_{{\rm GOA},\bot} (c,\omega ,b)\mathop{\cong}_{c\gg 1} (3\pi /4)R^{-1}(\omega )c^{-2}[{\rm erf}(2\pi b)/(2\pi b)]^{2}.\nonumber\\ \end{eqnarray} \tag{ 4.25 }$$

Assuming that b is fixed, the function F_GOA,⊥ depends on the pixel size h as follows, F_GOA,⊥ ∼ h². Therefore, for a fixed photon flux density in a projection image, N_ph,pxl/h², such that std_n,wht ∼ 1/h, the standard deviation of noise in a reconstructed slice is related to the pixel size as follows:

$$\begin{equation} {\rm std}_{n,\mu ,\bot} =\frac{{\rm const}}{\sqrt {N_{{\rm a}}} h}\frac{{\rm std}_{n,{\rm wht}}}{\sqrt {F_{{\rm GOA},\bot}}}\sim \frac{1}{h}. \end{equation} \tag{ 4.26 }$$

Regarding the corresponding gain function at large values of c, the following result is useful for simple quantitative assessment:

$$\begin{equation} G_{{\rm GOA},\bot} (c,\omega ,b)\mathop{\cong}_{c\gg 1} C_{b,\bot} R^{1/2}(\omega )c. \end{equation} \tag{ 4.27 }$$

where the proportionality constant C_b,⊥ depends only on the detector's parameter b, and in the case of a Gaussian MTF of the detector (with b = σ_det/(2h)) it is written as follows:

$$\begin{equation*} C_{b,\bot} =\left\{{\frac{1}{2(2\pi b)^{2}}-\frac{\exp [-(2\pi b)^{2}]}{\sqrt \pi {\rm erf}(2\pi b)(2\pi b)}} \right\}^{1/2}. \end{equation*}$$

For b = 0 the constant is C_0,⊥ = 0.577, while for b = 1/2 it is C_1/2,⊥ = 0.225.

In section 3.3, we discussed the choice of parameter (Δω/σ). It has been mentioned that at the x-ray energy of 20 keV (λ = 0.62 × 10⁻⁴ µm), the Si(1 1 1) double-crystal curve has the FWHM of 18.341 µrad. The corresponding value for σ is σ = FWHM/1.665 = 11.0 µrad. For the detector pixel size h of 1 µm and 10 µm, the corresponding values for (Δω/σ) are about 0.45 and 0.045, respectively.

Also we recall (see section 3.3 for details) that the ratio Δω/σ is energy independent (at least, at high x-ray energies, >20 keV).

Dependencies of the gain functions G_GOA,∥ and G_GOA,⊥ on the relative deviation of the analyser crystal, ω/σ, for various values of the parameter b, γ and Δω/σ, are shown in figure 7. It should be mentioned that, at similar imaging conditions (the same x-ray energy, detector and radiation dose), the 'parallel' imaging configuration results in a smaller gain compared to the 'perpendicular' configuration. This behaviour agrees with the above theoretical findings, namely that $G_{{\rm GOA},\vert \vert} \mathop{\sim}\limits_{c\gg 1} c^{1/2}$ while $G_{{\rm GOA},\bot} \mathop{\sim}\limits_{c\gg 1} \,c$. For example, in the case of γ = 10³, Δω/σ = 1/2 and b = 0, the maximum gain in the 'parallel' configuration is about 30, while in the 'perpendicular' configuration the maximum gain is more than 1000.

Zoom In Zoom Out Reset image size

On the other hand, the 'perpendicular' configuration is more sensitive to the detector's pixel size, h, through the angular spread Δω/σ, and to the ratio γ in the object. For instance, ten-fold reduction of Δω/σ, from 1/2 to 1/20, due to ten-fold increase of h, results in the reduction of the gain G by a factor of 10^1/2 and 10, for the 'parallel' and 'perpendicular' configuration, respectively. The resultant gain is about 9–10 for the 'parallel' configuration and is slightly bigger than 100 for the 'perpendicular' configuration, still an order of magnitude difference. The same results could be obtained if, instead of reducing the detector's pixel size, one would reduce the ratio γ by a factor of ten. Note that even in a less favourable case of γ = 10², Δω/σ = 1/20 and b = 0, the maximum gain in the 'parallel' configuration is about 3 while in the 'perpendicular' configuration the maximum gain is more than 10.

It should also be noted that the gain factor in both configurations is very sensitive to the width of the MTF of the detector, through the parameter b. Changing this parameter from 0 to 1/2 results in almost two-fold reduction of the gain factor in the 'parallel' configuration and in more than two-fold reduction of the gain factor in the 'perpendicular' configuration. Nevertheless, even in a much less favourable (but, probably, more realistic) case of γ = 10², Δω/σ = 1/20 and b = 1/2, the maximum gain in the 'parallel' configuration is about 1.7 while in the 'perpendicular' configuration the maximum gain is about 4.7.

The same numerical phantom and the imaging conditions as in section 3.4 have been used in numerical simulations. In total, 720 and 2880 projection images (in the angular range of 180°) have been calculated for the 2 cm and 8 cm phantoms, respectively. The resultant pixel size h in the projections was 10 µm. A single sinogram has been created for each projection data set. Poisson noise has been generated in the sinograms using a certain value for the mean absorbed dose delivered to the phantom. Two doses have been used in the numerical simulations, 5 and 500 mGy. The corresponding photon numbers per single detector pixel in an individual projection, in the beam incident on the object, as well as the average photon numbers per single pixel, in the beam incident on the detector, are summarized in tables 5 and 6.

For each sinogram, a processed copy was created using a 1D implementation of TIE-HOM phase-retrieval algorithm. All sinograms were subsequently processed using FBP CT reconstruction algorithm. The resultant reconstructed slices (2k × 2k and 8k × 8k pixels for 2 cm and 8 cm phantoms, respectively) contain β values and have been analysed by calculating the mean value and the standard deviation in the centre of the cylinder's cross-section and near its edge (using 256 × 256 and 1, 024 × 1, 024 regions for the 2 cm and 8 cm phantoms, respectively). The results of this analysis are presented in figures 8–11.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 8 shows measured values of the standard deviation of noise in the centre (red squares) and near the edge (blue circles) of the phantom, in the slices reconstructed using unprocessed contact projections. Figure 8 also shows the corresponding theoretical values for the standard deviation of noise (black solid line) calculated using equation (4.1) with ${\rm std}_{n,{\rm wht}}^{2} =N_{{\rm inc,pxl}}^{-1}$, where N_inc,pxl is the number of photons per single detector pixel (per single projection) in the beam incident on the object; these are given in tables 5 and 6. The theoretical values here completely ignore the attenuation of the incident beam by the object. As a result, they do not match the experimentally measured values. It is clear that in order to obtain accurate theoretical estimates for std_n,μ, proper correction for attenuation of the incident beam in the phantom is required. The following considerations have been used in order to derive a corresponding correction factor:

• (1)  
  the variance of noise in phase-contrast projection images is proportional to intensity, ${\rm std}_{n,z}^{2} \sim I_{z}$;
• (2)  
  phase contrast is usually small compared to intensity in the corresponding contact image, see equations (2.10),(2.22), so that in the noise analysis we can neglect phase contrast, ${\rm std}_{n,z}^{2} \sim I_{0}$;
• (3)  
  hence, the variance of noise in phase-retrieved projections is approximately proportional to the contact intensity, ${\rm std}_{n,0}^{2} \sim I_{0}$;
• (4)  
  sinograms contain not intensities but their negative logarithms, S₀ = −ln I₀, with the corresponding variance of noise written as follows, ${\rm std}_{n,S}^{2} \cong I_{0}^{-2} {\rm std}_{n,0}^{2} \sim I_{0}^{-1}$;
• (5)  
  backprojection step for each reconstructed pixel is essentially integration of the sinogram S along a sinusoidal trajectory, so that the variance ${\rm std}_{n,\mu}^{2}$ of noise in a pixel of the reconstructed slice can be approximated as follows (assuming 180° scan), ${\rm std}_{n,\mu}^{2} \sim F\equiv \pi^{-1}\int_0^\pi {{\rm d}\theta I_{0,\theta}^{-1} [r\sin (\theta -\theta_{0} )]}$,

where r and θ₀ are the polar coordinates of the pixel. In the case of a cylindrical object having the radius R and the linear attenuation coefficient μ, the correction factor, F_cyl, is calculated as follows:

$$\begin{equation*} F_{{\rm cyl}} =(2/\pi )\int_0^{\pi /2} {{\rm d}\theta \,\exp \{2\mu R[1-(r/R)^{2}\sin^{2}\theta ]^{1/2}\}} . \end{equation*}$$

This correction factor is a function of two dimensionless arguments, the ratio r/R and the product μR. In particular, for edge pixels (for which r/R ≈ 1) one has

$$\begin{eqnarray} &&\fl F_{{\rm cyl,e}} =(2/\pi )\int_0^{\pi /2} {{\rm d}\theta \,\exp (2\mu R\cos \theta )} \nonumber\\ &&=I_{0} (2\mu R)+L_{0} (2\mu R), \end{eqnarray} \tag{ 4.28 }$$

where I₀ is the modified Bessel function of the first kind and zeroth order and L₀ is the modified Struve function of the zeroth order. For central pixels (for which r/R ≈ 0) one has

$$\begin{equation} F_{{\rm cyl},c} =\exp (2\mu R). \end{equation} \tag{ 4.29 }$$

Correction factors, calculated using equations (4.28) and (4.29) for the two values of the diameter of the cylinder, are shown in table 7. Theoretical values of the standard deviation of noise, corrected for attenuation in the phantom using equations (4.28) and (4.29), are shown in figure 8 with red dashed and blue dashed–dotted lines, for central and edge pixels, respectively.

Analysis of figure 8 allowed us to make the following conclusions:

• (1)  
  As soon as the photon statistics in projection images is sufficient to approximate the Poisson probability distribution of noise with a Gaussian distribution, as, for example, in the case of 2 cm diameter and 500 mGy dose, a perfect match is observed between the experimentally measured and the theoretically predicted (using correction factors from equations (4.28) and (4.29)) values of std_β.
• (2)  
  In the opposite case, i.e. when the photon statistics in projection images is strongly Poissonian (the average number of photons per pixel is less than 10), as, for example, in the case of 8 cm diameter and 5 mGy dose (see tables 5 and 6 for details), no match is observed between the experimentally measured and the theoretically predicted values of std_β.

Figure 9 shows measured values of the standard deviation of noise in the centre (red squares) and near the edge (blue circles) of the phantom, in the slices reconstructed using phase-retrieved (using TIE-HOM) phase-contrast projections. Figure 9 also shows the corresponding theoretical values for the standard deviation of noise calculated using 1D version of equation (4.5) without correction for attenuation in the object (black solid line) and using correction factors (red dashed and blue dashed–dotted lines) given by equations (4.28)–(4.29). Analysis of figure 9 allows us to make the following conclusions:

• (1)  
  In the majority of considered cases (with only one exception, represented by the cylinder with 8 cm diameter and 5 mGy dose at 20 keV) noise statistics in the phase-retrieved projections is well approximated by the Gaussian distribution (as a consequence of weighted averaging over a large number of pixels during phase retrieval and the central-limit theorem). As a result, good match is observed between the experimentally measured and the theoretically predicted (and corrected for attenuation) values of std_β.
• (2)  
  Even when noise statistics in the phase-retrieved projections is not strictly Gaussian, as, for example, is the case for 8 cm diameter and 5 mGy dose at 20 keV, the relative error between the experimentally measured and the theoretically predicted values of std_β does not exceed 25%.
• (3)  
  Theoretical results presented in figure 9 match the corresponding experimental values much better than do the theoretical results presented in figure 8.

Figure 10 shows the average values of the imaginary part β of the refractive index of the cylinder, in the centre and at the edge of the reconstructed slices. Results shown in figure 10(a) correspond to reconstruction from conventional (contact) projections while results in figure 10(b) correspond to reconstruction from phase-retrieved (using TIE-HOM) propagation-based phase-contrast images. Analysis of figure 10 shows the following:

• (1)  
  PCT produces more accurate results than conventional CT.
• (2)  
  In the case of poor photon statistics, as for example in the case of the cylinder with 8 cm diameter and 5 mGy dose (for all energies in the range 20–45 keV), β values reconstructed using conventional CT are incorrect. This is due to the used sinogram creation procedure which assigned zeros to sinogram values, −lnI, for projection image pixels with zero intensity values. As a result, the sinograms and the β values, reconstructed using FBP, are strongly underestimated.
• (3)  
  PCT effectively solves the problem of poor photon statistics, stated in the previous item, since the phase-retrieval step results in local smearing of projection images and in significant reduction of the number of pixels with zero intensity. The only case in our numerical experiment that resulted in significant reconstruction errors was {8 cm, 5 mGy, 20 keV, central pixels} which had the poorest photon statistics.

Figure 11 shows the experimental (symbols) as well as the theoretical (solid black line) values for the gain factor, calculated using equation (4.10). Analysis of figure 11 indicates that the quality of match between the experimental and the theoretical gain values correlates with the quality of match between the experimental and the theoretical values for std_β, in conventional CT (see figure 8 and the related analysis). In particular, with the exception of the cases with very low photon statistics in projection images, the theoretically estimated gain factors are in a very good agreement with the corresponding experimentally measured gain factors. Importantly, there is no need for correction for attenuation in order to calculate the theoretical gain factors properly. This is because the correction factor in equation (4.23) equally affects the standard deviation in the slices reconstructed from conventional (contact) projection as well as the standard deviation in the slices reconstructed from the corresponding phase-retrieved phase-contrast projections.

We complete our numerical analysis by comparing the SNR values, defined as SNR = 〈β〉/std_n,β, in the PCT reconstructed slices calculated using the experimentally measured values of 〈β〉 and std_{n, β}, shown in figures 10 and 9, respectively, and the corresponding theoretical estimations. The SNR values in the two subregions of the reconstructed slices, i.e. at the centre of the slice and in a small area near the boundary of the cylinder, are compared in figure 12. Theoretical data points include the uncorrected values of SNR (black solid line), as well as the values corrected for the attenuation (red dashed and blue dashed–dotted lines).

Zoom In Zoom Out Reset image size

Analysis of figure 12 indicates a very good match of experimental and theoretical SNR values. The only exception is the case {8 cm, 5 mGy, 20 keV}, which corresponds to the poorest photon statistics in projection images. This misbehaviour correlates well with the results presented in figures 9 and 10. It should be emphasized that in the case of 8 cm diameter cylinder, without correction, theoretically calculated SNR values predict a wrong optimal energy (in the range 20–45 keV), namely, 20 keV. On the other hand, the experimental and the corrected theoretical results indicate that the optimal energy (out of six considered values) is 30 keV.

Having validated the theoretical formalism proposed in sections 4.2–4.4, we now present a quantitative comparison of the three imaging methods: conventional (absorption) imaging, PBI and ABI, the latter two are presented in conjunction with the corresponding phase-retrieval methods, TIE-HOM and GOA-HOM, respectively. These three methods are applied to the cylindrical phantom used throughout this paper. The comparison is carried out in terms of the SNR in the tomographic slices (containing the distribution of the linear attenuation coefficient μ) reconstructed using the FBP algorithm. The SNR is defined as follows, SNR = μ/std_n,μ, and the standard deviation of noise in the reconstructed slices is calculated using the analytical formulae derived in sections 4.1, 4.2 and 4.4. Figure 13 shows the SNR versus the x-ray energy, in the centre of tomographic slices of the phantom, reconstructed from TIE-HOM phase-retrieved propagation-based phase-contrast projections (b), (c) as well as from GOA-HOM phase-retrieved analyser-based phase-contrast projections (d),(e). In the case of PBI, either the propagation distance z was fixed to 100 cm (figure 13(b)) or the Fresnel number N_F = h²/(λ z) was fixed at 1 (figure 13(c)). In the case of ABI, results in 'parallel' and 'perpendicular' configurations are presented in figures 13(d) and (e), respectively. The detector pixel size was h = 10 µm and the mean absorbed dose was D_abs = 1 mGy.

Zoom In Zoom Out Reset image size

It should be noted that in contrast to the results presented in figure 12 (where 1D phase retrieval has been carried out), figure 13(b) shows SNR in the case of 2D phase retrieval. As a result, even with smaller dose, 1 mGy (figure 13(b)) instead of 5 mGy (figures 12(b) and (d)), one obtains better SNR in the reconstructed slices. Indeed, according to table 8, approximately a 3-fold increase of the gain factor is observed in the case of the 2D phase retrieval, compared to the 1D phase retrieval. This results in the same improvement in SNR, which increases by approximately a factor of 3. On the other hand, a 5-fold dose reduction results in 5^1/2-fold reduction in the SNR. Thus approximately a 1.3-fold overall improvement of SNR is observed in figure 13(b) in comparison to figures 12(b) and (d).

Analysis of figure 13 shows that, regardless of the imaging modality and imaging conditions, the optimum energy, that maximizes SNR, increases with the object's size. At the same time, the corresponding maximal SNR decreases with the object's size. Note, however, that SNR reduction in the case of CT reconstruction is very modest in comparison with the corresponding SNR reduction observed in individual projection images (see figure 4 for details).

Under the considered conditions, namely, the free-space propagation distance z = 100 cm and a perfect detector (b = 0) with the pixel size h = 10 µm, the phase-contrast CT reconstruction using TIE-HOM phase retrieval outperforms that with GOA-HOM phase retrieval. Moreover, the popular 'parallel' configuration in ABI produces approximately an order of magnitude worse results compared to PBI CT, in terms of SNR, as follows from figure 13, as well as from table 8 (see gain values there).

While TIE-HOM (with fixed propagation distance) and GOA-HOM in 'perpendicular' configuration show similar performances, in terms of SNR, it is worth mentioning that GOA-HOM requires higher x-ray energies for optimal performance, about 10–15 keV higher compared to that for TIE-HOM. For instance, in the case of the cylinder with the diameter of 8 cm, the optimum energies in TIE-HOM and GOA-HOM case are 30 keV and about 40 keV, respectively. The use of high x-ray energies in ABI often involves significant technical complications; in particular, issues with the x-ray beam stability and the flux.

In order to improve the SNR in reconstructed slices one could consider increasing the detector pixel size, h. Note, however that according to our findings in sections 4.2 and 4.5, in the case of PBI with subsequent phase retrieval using TIE-HOM, the SNR is almost independent of the pixel size (as long as the free-space propagation distance is fixed), while the corresponding gain factor is approximately inversely proportional to the second power of h. In the case of ABI (in 'perpendicular' configuration) with subsequent phase retrieval using GOA-HOM, the SNR is proportional to the pixel size while the corresponding gain factor is inversely proportional to it. This means, in particular, that (i) ABI(⊥) + GOA-HOM can outperform PBI(z) + TIE-HOM for sufficiently large detector pixels; (ii) however, increasing the detector pixel size, one decreases the corresponding gain factor compared to conventional contact imaging and, eventually, there may be no benefit in using phase-contrast projections in place of conventional (contact) projection images for large detector pixel sizes.

Figure 13(c) shows the SNR in the slices reconstructed using TIE-HOM, but instead of fixing the free-space propagation distance, as in figure 13(b), the Fresnel number N_F = h²/(λz) = 1 is fixed, thus providing maximum performance of the TIE-HOM phase retrieval for all x-ray energies. This approach results in a significant increase of the optimum x-ray energy and a noticeable improvement in the SNR, compared to the case with a fixed propagation distance shown in figure 13(b). Moreover, the optimum energy here is similar to that in the case of ABI with GOA-HOM phase retrieval in 'perpendicular' configuration shown in figure 13(e). Also, when the Fresnel number is fixed, the corresponding SNR in the slices reconstructed using TIE-HOM phase-retrieved projections is proportional to the second power of the pixel size (exactly as in the case of CT reconstruction from contact projection images). This indicates that by increasing the pixel size, the difference between the maximum SNR in the slices reconstructed using PBI + TIE-HOM(N_F = const) and that in the slices reconstructed using ABI + GOA-HOM(⊥) can only increase. On the other hand, one might expect a significant improvement of the relative performance of ABI CT (with respect to PBI CT) in the case of small detector pixels.

In the case of PBI there is another possibility to improve SNR in the reconstructed slices, by increasing the free-space propagation distance, z, and using CTF-HOM phase-retrieval algorithm instead of TIE-HOM. We recall, however, that there is a limit for the corresponding gain, which was estimated to be around 0.3γ (see section 4.3 and figure 6, in particular). Applying this result to the case of the cylindrical phantom with the diameter of 8 cm at the optimum energy of 30 keV (see figure 13(b)), we obtain that the ratio γ for the breast tissue is 2332 and hence the upper limit on the gain in propagation-based phase-contrast CT is about 700, regardless of the pixel size (!). Note, however, that this might involve very large propagation (object-to-detector) distances, z ∼ h²/λ. For instance, at 30 keV (λ ≈ 0.4 × 10⁻⁴ µm) and h = 10 µm, the corresponding distance is z ≈ 2.5 m and is relatively easy to manage. For a larger pixel size, h = 20 µm, the corresponding distance becomes approximately 10 m, which can only be realized at specialized (large) imaging facilities, such as synchrotrons.

Also, we would like to emphasize the fundamental difference, in terms of the relationship between the variance of noise and the photon statistics, of the three imaging methods, including the conventional radiography, PBI combined with TIE-HOM phase retrieval and ABI combined with GOA-HOM phase retrieval. These relationships follow from our theoretical analysis presented in sections 4.2 and 4.4 as well as in this section, and are summarized in table 9. Also shown in table 9 are the typical gain values for the two PCT methods.