1. Introduction

Due to the advent and accessibility of synchrotron and laboratory X-ray sources [1, 2, 3] that possess high degrees of coherence, new imaging methods, broadly characterized as X-ray phase-contrast imaging methods [4, 5, 6, 7, 8, 9, 10, 11, 12, 13] are being developed that have dramatic advantages over conventional radiographic and tomographic X-ray imaging systems. Unlike conventional radiographic methods, phase-contrast imaging methods exploit a contrast mechanism based on differences in the complex-valued X-ray refractive index values of tissue. At diagnostic X-ray energies, variations in the real component of the refractive index of tissues are several orders of magnitude larger than are variations in the imaginary component, or equivalently, the X-ray attenuation coefficient[14]. Consequently, phase-contrast imaging may permit the visualization and quantification of tissues that have identical, or very similar, X-ray absorption properties. Additionally, unlike absorption-contrast, the phase-contrast mechanism persists at relatively high energies, which permits low-dose imaging[15]. Experimental studies of X-ray phase-contrast imaging have revealed significantly enhanced contrast and tissue discriminability over conventional radiographic methods in applications including cancer detection [16, 17, 18, 19] and cartilage imaging [20, 21].

Most experimental implementations of X-ray phase-contrast imaging [4] that hold promise for clinical imaging are broadly characterized as being either analyzer-based or propagation-based. In analyzer-based methods [22, 11], the object is irradiated with a quasi-monochromatic, parallel-beam, X-ray wavefield, and an analyzer crystal or other X-ray diffracting element [13] is present between the object and detector. The analyzer crystal diffracts the components of the transmitted wavefield traveling at or near the crystal’s Bragg angle, thereby rejecting all wavefield components traveling outside a narrow angular range. This angular filtering of the transmitted wavefield produces a phase-contrast enhancement in the recorded intensity image. Propagation-based methods [23, 6, 24, 25], also known as in-line methods, do not require the use of diffractive X-ray elements, and employ a classic in-line holographic measurement geometry that is essentially similar to magnification radiography [9]. The object is irradiated with a plane-wave or paraxial X-ray wavefield possessing a sufficient degree of spatial coherence, and the intensity of the transmitted wavefield is recorded by a detector placed at some non-zero distance behind the object. In this case, the measured intensity distribution represents an in-line (Gabor) hologram [26] that contains coded information about the refractive properties of the object. The phase shifts that are introduced into the probing wavefield by the refractive index variations within the object are transferred into intensity variations in the final measurement by the process of free-space (Fresnel) wave propagation between the object and detector.

Quantitative X-ray phase-contrast imaging methods seek to reconstruct separate images that depict the object’s projected absorption and real-valued refractive index distributions, which reflect two distinct and complementary intrinsic object properties. To achieve this, the measured intensity data must generally be recorded at two or more distinct “states” of the imaging system [27]. For example, in analyzer-based methods, different measurement states correspond to distinct orientations of the analyzer crystal, while in propagation-based methods they could correspond to distinct object-to-detector distances. Quantitative phase-contrast imaging methods are computed-imaging methods, and require use of reconstruction algorithms for image formation. When implemented in tomographic mode, an estimate of the three-dimensional (3D) complex refractive index distribution can be reconstructed by use of phase-contrast tomography methods [28, 29, 30, 31, 21].

Under certain assumptions regarding the transmitted wavefield, phase-contrast imaging systems can be modeled as linear, shift-invariant, imaging systems. Assuming measurements are acquired at two distinct states of the imaging system, Fourier-based reconstruction formulas can be derived [27, 32, 33, 34] from knowledge of the imaging system’s optical transfer function. A common feature of these reconstruction formulas, however, is the presence of isolated Fourier domain singularities. In practice, the Fourier components of the projected object properties residing near the singularities can contain greatly amplified noise levels, resulting in noisy and distorted images. Because the locations of the singularities are determined by the chosen two measurement states, their effects can be mitigated by acquiring intensity measurements at multiple (>2) measurement states [28, 35]. For example, different pairs of intensity measurements can be utilized to reconstruct different Fourier components of the object properties, which provides the opportunity to potentially avoid the singularities. However, such simple strategies do not exploit fully the statistically complementary information contained in the available measurement data, and can result in images with suboptimal statistical properties.

In this work, we propose and investigate a methodology for image reconstruction in quantitative phase-contrast imaging when multiple (>2) intensity measurements are available. Linear estimators are proposed that combine the available intensity measurements to produce statistical estimates of the projected object properties whose Fourier components possess optimally reduced variances. This general strategy is inspired by a recent study of intensity diffraction tomography [36] by our group. Explicit forms of the reconstruction formulas are derived for propagation-based phase-contrast imaging, where a general noise model and finite-sampling effects are considered. Computer-simulation studies are conducted to demonstrate the efficacy of the proposed method.

2. Background

In this section, the salient features of X-ray phase-contrast imaging are reviewed. We refer the reader to the monograph by Paganin [37] for a comprehensive treatment of phase-contrast image formation.

2.1. Interaction of X-ray wavefield with object

As depicted in Fig. 1, consider that an object is irradiated by a monochromatic X-ray wavefield U_i with wavelength λ, which is traveling along the positive z-axis. The effects of imperfect wavefield coherence will not be considered, but can be addressed as in [38, 39]. The object is characterized by its complex-valued refractive index distribution

where $\overrightarrow{r}$ = (x,y,z), and δ($\overrightarrow{r}$) and β($\overrightarrow{r}$) quantify the X-ray refractive and absorption properties of the object. The quantity β($\overrightarrow{r}$) is related to the linear X-ray attenuation coefficient µ($\overrightarrow{r}$) as

where $k\equiv \frac{2\pi }{\lambda }$ is the wavenumber. Note that classic radiographic methods are sensitive only to variations in β($\overrightarrow{r}$), while phase-contrast methods are sensitive to variations in both δ($\overrightarrow{r}$) and β($\overrightarrow{r}$). The wavefield U_o(x,y) on the object plane, which has been transmitted through the object, is given by

where T(x,y) is the transmission function that can be expressed generally as

The amplitude modulus M(x,y)=exp[-A(x,y)] and the phase shift ϕ(x,y) describe how the amplitude and wavefront (i.e., phase), respectively, of the probing wavefield are perturbed by the presence of the object. They are related to imaginary and real components, respectively, of n($\overrightarrow{r}$) as

where the line-integrals are computed over the support of the object.

 

Fig. 1. A schematic of a generic X-ray phase-contrast imaging system.

Download Full Size | PDF

2.2. Linear shift-invariant X-ray phase-contrast imaging systems

Consider again the generic X-ray phase-contrast imaging system shown in Fig. 1. Let U_m(x,y) denote the transmitted wavefield on a detector plane of constant z, which is downstream from the object plane. The integer-valued subscript m is employed to denote the state of the imaging system. For many analyzer- and propagation-based phase-contrast imaging systems, U_m(x,y) and U_o(x,y) can be regarded as the output and input, respectively, of a linear and shift-invariant coherent imaging system [40, 41]. In this case,

where ∗ denotes a two-dimensional (2D) convolution and G_m(x,y) describes the impulse response of the system in state m. In propagation-based imaging, for example, G_m(x,y) describes the Fresnel propagator [26], with distinct values of m corresponding to different object-to-detector distances. Alternatively, in analyzer-based imaging, G_m(x,y) describes the coherent impulse response of the analyzer crystal or other diffractive element(s) employed, with distinct values of m corresponding to different orientations of the analyzer. In hybrid systems [39], G_m(x,y) describes the net effect of both.

On the detector plane, the intensity I_m(x,y)=|U_m(x,y)|² is recorded, which represents a radiograph with mixed absorption- and phase-contrast. From knowledge of the measured intensity I_m(x,y), a modified data function can be defined as

where I_i=|U_i|² is the intensity of the incident X-ray beam. Let K̃_m(u,v) denote the 2D Fourier transform (FT) of K_m(x,y) defined as

In biomedical imaging applications, the conditions of |A(x,y)|≪1 and slowly varying ϕ(x,y) can often be met [10, 42]. Under such conditions, it can be shown that [43]

where G̃^a_m(u,v) is the amplitude transfer function (ATF):

and G̃^p_m(u,v) is the phase transfer function (PTF):

Here, G̃_m(u,v), Ã(u,v), and ϕ̃(u,v) are the 2D FTs of G_m(x,y), A(x,y), and ϕ(x,y), respectively, and G̃*_m(·, ·) denotes the complex conjugate of G̃_m(·, ·). The interested reader is referred to Ref. [43] for a detailed derivation of the imaging model in Eq. (9).

Equation (9) relates the measured intensity I_m(x,y), or equivalently K_m(x,y), to the 2D FTs of the sought-after quantities A(x,y) and ϕ(x,y). If an additional measurement I_n(x,y) is obtained when the imaging system is in state n≠m, Ã(u,v) and ϕ̃(u,v) can be determined algebraically as

Equations (12a) and (b) represent reconstruction formulas for quantitative phase-contrast imaging. From the determined Ã(u,v) and ϕ̃(u,v), A(x,y) and ϕ(x,y) are computed by application of the inverse 2D FT. Note that Eqs. (12a) and (b) contain isolated poles at frequency components (u,v) for which

When Eq. (12) is applied to noisy, or otherwise inconsistent, measurement data, the Fourier components (u,v) of Ã(u,v) and ϕ̃(u,v) residing near the poles will contain greatly amplified noise levels. This can result in noisy and/or distorted images. In the remainder of this article, we describe a statistically principled method for circumventing this when measurements corresponding to multiple (>2) states of the system are available.

3. Variance reduction in quantitative X-ray phase-contrast imaging

We consider that intensity measurements I_m(x,y) are acquired at three distinct states m=1,2,3 of the system. The results below are generalized to the case of an arbitrary number of measurements in the Appendix. The intensity data function I_m(x,y) is interpreted as a stochastic process, which reflects that the measurements are contaminated by stochastic errors such as detector noise. Our goal is to exploit the statistically complementary information in the available measurements to reduce the variance of the estimated Ã(u,v) and ϕ̃(u,v), and thereby mitigate the large amplification of noise due to poles in the reconstruction formulas.

Because the reconstruction formulas in Eq. (12) require intensity measurements at two distinct states, $\left(\begin{array}{c}N\\ 2\end{array}\right)$ estimates of Ã(u,v) and ϕ̃(u,v) can be computed from knowledge of measurements obtained at N system states. When reconstructed from noisy measurements, these estimates will be generally distinct. The notation à _m,n(u,v) and ϕ̃_m,n(u,v), m≠n=1,2,3, will be employed to describe the estimates for the case N=3, where the subscripts denote that measurements I_m(x,y) and I_n(x,y) were employed. Because the locations of poles in Eq. (12) depend on the choice of measurement states, the components of à _m,n(u,v) or ϕ̃_m,n(u,v) that are highly contaminated by noise will be determined by the measurement state pair (m,n).

A natural strategy for mitigating noise amplification is to combine the à _m,n(u,v) or ϕ̃_m,n(u,v) in a way that attempts to cancel the poles in each two-state estimate [36], to form final estimates Ã(u,v) or ϕ̃(u,v), respectively, that possess reduced variances for all (u,v). However, it should be noted that only two of the three available estimates are independent in the sense that

with l≠m≠n=1,2,3. In other words, any estimate can be expressed as a linear combination of the remaining two. The coefficients in Eqs. (14) and (15) are frequency-dependent and given by

It can be verified that

Accordingly, final estimates of Ã(u,v) and ϕ̃(u,v) that exploit statistically complementary information in the three intensity measurements can be formed as

where the combination coefficients satisfy

Equations (18) and (19) represent linear estimators for producing unbiased estimates of ϕ̃(u,v) and Ã(u,v), respectively. Because the combination coefficients ω^ϕ _m,n(u,v) and ω^a _m,n(u,v) are frequency-dependent, they can be designed to cancel poles present in the ϕ̃_m,n(u,v) and à _m,n(u,v), respectively. Moreover, as described next, they can be designed to optimally reduce the variance of ϕ̃(u,v) and Ã(u,v). We consider first the problem of producing estimates ϕ̃(u,v) having reduced variances. The following notation will be employed:

where ‘Var’ and ‘Cov’ denote the variance and covariance, respectively of a random process. The variance of ϕ̃(u,v) is given by

where the superscript * denotes the complex conjugate. In order to minimize Var{ϕ̃(u,v)}, we need that

where R _1,2 and I _1,2 are the real and imaginary components of ω^ϕ _1,2, respectively. The solution of these equations yields

The optimal choice of ω^ϕ _1,3(u,v) is subsequently determined by Eq. (20a).

The specification of the combination coefficients ω^a _1,2(u,v) and ω^a _1,3(u,v) that optimally reduce the variance of Ã(u,v) via Eq. (19) are also determined by Eq. (26) when the quantities in Eqs. (21)–(23) are redefined appropriately. A heuristic method for choosing the combination coefficients that can effectively mitigate noise amplification when the noise model is not known is described later.

As described by Eq. 12, the reconstruction formulas for determining the phase function ϕ(x,y) and attenuation function A(x,y) are described by simple algebraic forms in the Fourier domain. Consequently, the large amplification of noise due to poles in the reconstruction formulas can be mitigated in a mathematically straightforward and physically understandable way in the Fourier domain. Reducing the Fourier domain variance of the phase and attenuation estimates generally leads to spatial domain estimates that possess reduced variances. This can be understood by noting that ${\iint }_{{\mathbf{R}}^2}\text{d}x\text{d}y\mathrm{Var}\left\{\varphi (x,y)\right\}={\iint }_{{\mathbf{R}}^2}dud\nu \mathrm{Var}\left\{\tilde{\varphi }(u,\nu )\right\}$ The left hand-side of this equation defines the global variance of the sought-after phase function. This indicates that a phase function possessing a reduced global variance can be obtained from an estimate ϕ̃(u,v) with a reduced variance. Because Var{ϕ(x,y)} is nonnegative, a lower global variance suggests, in general, lower local variances in the determined phase function. The same observation holds true for A(x,y).

4. Application to multi-plane propagation-based imaging

In the remainder of the article, the general methodology described in Section 3 is investigated within the context of propagation-based X-ray phase-contrast imaging. In this section, the explicit forms of the optimal estimators in Eqs. (18) and (19) are derived in their continuous forms. The effects of finite sampling are examined in Section 5.

4.1. Two-state reconstruction formulas

In propagation-based imaging, the different system states can correspond to different objectto- detector distances. We consider that three intensity measurements I_m(x,y) are acquired on distinct detector planes z=z_m, where m=1,2,3. The impulse response in Eq. (6) corresponds to the Fresnel propagator

and its 2D FT yields the transfer function

The corresponding reconstruction formulas [24] are found by use of Eq. (12):

where f ²≡u ²+v ² and

As before, the indices m, n satisfy m=1,2, n=2,3 with n>m.

Note that D _m,n=0 is equivalent to Eq. (13), and specifies the locations of poles in the reconstruction formulas. Equation (29) contains poles corresponding to frequencies (u,v) that satisfy

where l is an integer. One such pole is located at zero-frequency u=v=0, indicating that the low-frequency components of ϕ̃(u,v) will contain highly amplified noise levels. The existence of additional poles, away from the origin of the Fourier space, depends on the detector resolution and the detector pair spacing. Let (u_M,v_M) denote the maximum spatial frequencies recorded by the detector. Additional poles in the reconstruction formulas described by Eqs. (29) and (30) will be present when

Equation (33) indicates that for a fixed detector resolution, additional poles in ϕ̃_m,n(u,v) will emerge when the detector spacing (z_m-z_n) is sufficiently large. Likewise, this discussion of poles also applies to à _m,n(u,v), with the exception that the pole at u=v=0 is not present due to a cancellation.

4.2. Second-order statistics for determination of optimal combination coefficients

In order to determine the optimal combination coefficients ω^ϕ _1,2(u,v) and ω^ϕ _1,3(u,v) that minimize the variance ϕ̃(u,v) [via Eq. (18)], the variance and covariance information in Eqs. (21) and (22) must be determined. Knowledge of the analogous quantities involving Ã(u,v) is required for determination of the optimal combination coefficients of ω^a _1,2(u,v) and ω^a _1,3(u,v) that minimize the variance of Ã(u,v).

From Eqs. (29) and (30), one finds readily that

The joint covariances of the Fourier data are computed as

It should be noted that Eqs. (36) and (37) are real-valued, and therefore the imaginary components of the optimal combination coefficients will vanish [see Eq. (26b)].

4.3. Heuristic determination of combination coefficients

When the second-order statistics in Section 4.2 are not known, the optimal combinations coefficients ω^a _m,n(u,v) and ω^ϕ _m,n(u,v) cannot be computed. However, the large noise amplification due to poles in ϕ̃_m,n(u,v) and à _m,n(u,v) can still be effectively mitigated by suitable heuristic specification of the combination coefficients. From Eqs. (34)–(35), near the locations of poles, we find that

Heuristic combination coefficients ω ^heur _m,n (u,v) should have a (u,v)-dependence that is inversely proportional to that indicated in Eqs. (38)–(39). Moreover, at locations of poles we should have ω ^heur _m,n (u,v)≡0. Based on these requirements, the ω ^heur _m,n (u,v) for use in estimating ϕ̃(u,v) via Eq. (18) can be chosen as

where α^ϕ _1,2 and α^ϕ _1,3 are defined in Eq. (16). When estimating Ã(u,v), α^ϕ _1,2 and α^ϕ _1,3 are replaced by α^a _1,2 and α^a _1,3, respectively. Due to their construction, ω ^heur _1,2 (u,v) and ω ^heur _1,3(u,v) satisfy the normalization condition in Eq. (20).

5. Computation of optimal combination coefficients with consideration of finite sampling

Below, the second-order statistics described in Section 4.2, and subsequently the optimal combination coefficients ω^a _m,n(u,v) and ω^ϕ _m,n(u,v) are computed explicitly with consideration of finite sampling effects.

5.1. Noise model and finite sampling considerations

The discretely sampled intensity data on a measurement plane z=z_m is denoted as

where r and s are integer-valued indices that reference detector elements, and Δx=Δy denotes the element dimension in a square detector array of dimension L×L. Equation (41) assumes idealized (Dirac delta) sampling, namely, the averaging effects of sampling aperture are not considered. However, the analysis follows can be generalized to address such effects. The square bracket, ‘[·]’, is introduced to represent the functions whose arguments are discretely sampled. We assume the noise model satisfies [35, 44]

where I ⁰ _m[r,s] denotes the noiseless intensity data and the signal-dependent noise n_m(r,s) has zero-mean and

where the real-valued quantity σ ²(z_m) can depend on the detector location z_m [35]. We also assume

where δ_mn denotes the Kronecker delta function.

5.2. Second-order statistics of discrete data functions

In order to compute the second-order statistics described in Section 4.2, knowledge of the variance of Ĩ(u,v) is required. To compute this from knowledge of discretely sampled data, the continuous 2D FT will be approximated by use of the discrete Fourier transform (DFT) [45].

The 2D DFT [45] of Eq. (42) can be written as follows

where

with p, q denoting the integer-valued Fourier indices that are conjugate to [r,s], and N specifying the number of detector elements in each dimension of the square 2D detector. The variance of ñ_m[p,q] is computed as

where ‘E{·}’ denotes the statistical expectation operator. The continuous and discrete FTs of the intensity data are related as [45]

where $\Delta u=\Delta v=\frac{1}{L}$ are the frequency domain sampling intervals along u and v axes. By use of Eqs. (49) and (50), we find that

Finally, using Eq. (51) with the results of Section 4.2, we arrive at

where, as before, m=1,2, n=2,3 with n>m, and

5.3. Explicit forms of optimal combination coefficients

By use of Eqs. (52)–(55) and Eq. (26), we arrive at the optimal combination coefficients for use in estimating ϕ̃(u,v) or Ã(u,v) via Eq. (18) or (19), respectively:

where

and

where σ ² _m≡σ ²(z_m). The corresponding forms of ω^ϕ _1,3(u,v) and ω^a _1,3(u,v) are determined by Eq. (20).

If ∑^N-1 _r=0 ∑^N-1 _s=0(I ⁰[r,s;z_m])² does not vary significantly as a function of m (i.e., the detector location z_m), Eq. (56) can be expressed in the somewhat simplified form:

6. Numerical Studies

A preliminary computer-simulation study of propagation-based X-ray phase-contrast imaging was conducted to corroborate and quantitatively investigate the proposed reconstruction methods.

6.1. Numerical phantom and in-line measurement geometry

The in-line measurement geometry shown in Fig. 2 was assumed. A monochromatic X-ray plane-wave with wavelength λ=1×10^-10 m, propagated along the z-axis and irradiated an object. Three detector planes located at z=z_m, m=1,2,3, were considered to be behind the object. The detector contained 1024×1024 elements of dimension of 1 µm², and was assumed to have otherwise idealized physical properties. Two measurement geometries were considered, which will be referred to as Geometry ‘A’ and Geometry ‘B’. In Geometry ‘A’, the detector planes were positioned at z ₁=19 mm, z ₂=96 mm, z ₃=182 mm, while the corresponding positions in Geometry ‘B’ were z ₁=12 mm, z ₂=38 mm, and z ₃=72 mm.

 

Fig. 2. The measurement geometry of propagation-based X-ray phase-contrast imaging employing multiple detector-planes.

Download Full Size | PDF

A 3D mathematical phantom comprised of five uniform ellipsoids possessing different complex-valued X-ray refractive index values representative of soft tissue was employed to represent the object to-be-imaged. The semi-axes of the largest ellipsoid were 188.416 µm, 163.840 µm, and 141.312 µm. From knowledge of the phantom, the projected object properties ϕ(x,y) and A(x,y) were computed according to Eq. (5) and are displayed in Fig. 3.

 

Fig. 3. Images of the true object properties (a) A(x,y) and (b) ϕ(x,y).

Download Full Size | PDF

6.2. Measurement data and simulation studies

The simulated intensity measurements were computed as follows. From knowledge of ϕ(x,y) and A(x,y), the transmitted wavefield U_o(x,y) on the object plane was computed according to Eqs. (3) and (4). Subsequently, sampled values of the wavefieldU_m(x,y) on each detector plane z=z_m, m=1,2,3, was computed by use of Eq. (6) with G_m(x,y) specified by the Fresnel propagator in Eq. (27). The convolution in Eq. (6) was computed by use of the 2D fast Fourier transform (FFT) algorithm. The intensity data I_m[r, s] were then computed as the square of the wavefield modulus on each detector plane.

Noisy measurement data I_m[r,s] were computed according to Eq. (42). The noise process n_m[r,s] was described by a uncorrelated Gaussian distribution whose variance was determined by Eq. (43) with σ ²(z_m)=0.05%. Ensembles of 1000 noisy realizations of I_m[r,s] were computed for m=1,2,3.

6.3. Image reconstruction

Estimates ϕ̃_m,n(u,v) and à _m,n(u,v) were computed from each pair of noisy intensity data by use of Eqs. (29) and (30), respectively. The presence of poles can pose considerable difficulty in determining these estimates. Simply setting ϕ̃_m,n(u,v)=0 or à _m,n(u,v)=0 at the locations of poles will lead to inaccuracy in the resulted estimates. Additionally, even if the poles are avoided, the data errors can be greatly amplified in the vicinity of poles, where the denominators of Eqs. (29) and (30) take on small values. In current studies, the reconstruction formulas were regularized by setting the estimates of ϕ̃_m,n(u,v) and à _m,n(u,v) to zeros in the vicinity of poles when D _m,n≤2×10^-7. These estimates were combined, according to Eqs. (18) and (19), to form final estimates ϕ̃(u,v) and Ã(u,v) that possess optimally reduced variances. The required combination coefficients were computed according to Eqs. (56) and (57). Because σ ²(z_m) was fixed at a constant value, it can be verified that, in this special case, the optimal combination coefficients given in Eqs. (56) and (57) are identical to the heuristic ones defined by Eq. (40). Corresponding estimates of ϕ(x,y) and A(x,y) were computed by application of the 2D inverse FFT algorithm. The variances of the reconstructed object properties in both the Fourier and spatial domains were estimated empirically.

6.4. Numerical results

Numerical studies were conducted to corroborate the noise analysis described in Section 5. Specifically, the theoretically predicted Var{ϕ̃(u,v)}, as stated in Eq. (24), was compared to an empirically determined estimate. The same was done for Var{Ã(u,v)}. When computing the theoretically predicted Var{ϕ̃(u,v)}, Eqs. (34), (36), and (56) were employed with Eq. (24). Similarly, Eqs. (35), (37), and (57) were employed to determine the theoretically predicted Var{Ã(u,v)}. The empirical estimate of Var{ϕ̃(u,v)} was determined as follows. Firstly, empirical estimates of the two-state variance and covariance functions in Eqs. (21) and (22) were computed from the ensembles of noisy intensity data. Secondly, these quantities were employed to determine estimates of the optimal combination coefficients ωϕ_m,n(u,v) via Eq. (26). Lastly, an empirical estimate of Var{ϕ̃(u,v)} was determined from an ensemble of noisy images reconstructed by use of Eq. (18). The same procedure was followed for determining the empirical estimates of Var{Ã(u,v)}. Figures 4 and 6 display the determined variance maps corresponding to Geometry ‘A’ and Geometry ‘B’, respectively, which have been logarithmically transformed for display purposes. In each figure, subfigures (a) and (b) display the theoretically predicted and empirically determined images of log [Var{Ã(u,v)}], respectively, while the corresponding images of log [Var{ϕ̃(u,v)}] are contained in subfigures (c) and (d), respectively. The theoretically predicted and empirical variance maps appear nearly identical. This is confirmed by Figs. 5 and 7, in which horizontal profiles through the centers of the theoretically predicted and empirical variance maps are superimposed, respectively. These results demonstrate excellent agreement for all Fourier components.

Figure 8 displays estimates of A _m,n(x,y) reconstructed from noisy intensity data measured in Geometry ‘A’ by use of detector planes (a) (1,2), (b) (1,3), (c) (2,3). The ‘optimal’ estimate of A(x,y) obtained by use of Eq. (19), which employs all three intensity measurements, is shown in (d). The corresponding estimates of ϕ _m,n(x,y) and ϕ(x,y) are shown in subfigures (e)–(g) and (h), respectively. The excessively noisy appearances of the A _m,n(x,y) in subfigures (a)–(c) are due to the low absorption contrast of the object. In this measurement geometry, the relatively large detector spacings result in the occurrence of extra poles away from the origin in Fourier space [see Eq. (33)]. This creates severe noise amplification in the high-frequency components of the two-state reconstructions A _m,n(x,y) and ϕ _m,n(x,y). The optimal estimates of A(x,y) and ϕ(x,y) shown in subfigures (d) and (h), respectively, contain obviously reduced noise levels as compared to the two-state reconstructions. The estimates of ϕ _m,n(x,y) and ϕ(x,y) appear to be contaminated by low-frequency noise, as evident by their lumpy background appearances. This is explained by the fact that all three estimates of ϕ _m,n(x,y) have a pole at the origin of Fourier space and therefore this pole cannot be removed by the estimator in Eq. (18). Note that taking a simple average of the available two-state estimates of ϕ _m,n(x,y) or A _m,n(x,y) would not be an effective reconstruction strategy because it does not mitigate large noise amplifications due to poles in the reconstruction formulas away from the origin of Fourier space.

 

Fig. 4. Images of theoretical and empirical estimates of Var{Ã(u, v)} measured in Geometry ‘A’ are displayed logarithmically in subfigures (a)–(b), respectively. The corresponding variance maps of ϕ̃(u, v) are contained in subfigures (c)–(d), respectively.

Download Full Size | PDF

 

Fig. 5. Variance profiles of images in Fig. 4. Subfigure (a) contains the theoretically and empirically determined variance profiles of Ã(u, v), which are depicted by solid and dashed curves, respectively. The corresponding variance profiles of ϕ̃(u, v) are shown in subfigure (b).

Download Full Size | PDF

 

Fig. 6. Images of theoretical and empirical estimates of Var{Ã(u, v)} measured in Geometry ‘B’ are displayed logarithmically in subfigures (a)–(b), respectively. The corresponding variance maps of ϕ̃(u, v) are contained in subfigures (c)–(d), respectively.

Download Full Size | PDF

 

Fig. 7. Variance profiles of images in Fig. 6. Subfigure (a) contains the theoretically and empirically determined variance profiles of Ã(u, v), which are depicted by solid and dashed curves, respectively. The corresponding variance profiles of ϕ̃(u, v) are shown in subfigure (b).

Download Full Size | PDF

Variance data for the images contained in Fig. 8 are shown in Fig. 9. Figures 9(a) and (b) display empirical variance estimates of the reconstructed Fourier components. Each subfigure contains profiles corresponding to the Fourier variance estimated by use of detector pair (1,2) (dashed curve), detector pair (1,3) (dashdotted curve), detector (2,3) (dotted curve), and the optimal ones reconstructed by use of Eqs. (18) or (19) (solid curve). As expected, the optimal estimates possess variances that are lower than any of the two-state estimates, for all frequency components. Subfigures (c) and (d) display empirical estimates of the corresponding images in the spatial domain. The labeling of the profiles is the same as described above. The variances of the optimally estimated images are seen to be lower than the two-state reconstructions, for all values of (x,y). This reflects that estimates having reduced Fourier variances will generally have reduced variances in the spatial domain.

The corresponding results obtained with Geometry ‘B’ are shown in Figs. 10 and 11. From Fig. 10, we find again that the optimal estimates of A(x,y) and ϕ(x,y) contain obviously reduced noise levels as compared to the two-state reconstructions. The quantitative variance data in Fig. 11 confirms that the optimal estimates of A(x,y) and ϕ(x,y) contain lower variances than the two-state reconstructions.

In practice, the detector locations cannot be determined with perfect precision. The uncertainty in the detector positions will introduce inconsistencies between the measured data and assumed imaging model. In order to investigate the effects of this, images were reconstructed from the noisy intensity measurements corresponding to Geometry ‘A’, where imperfect knowledge of the detector positions was assumed. The assumed detector positions ẑ_m were related to the true locations zm by ẑ_m=z_m+ε(z_m), with m=1,2,3, where ε(z_m) represents the positioning errors. The three sets of errors ε(z_m) contained in Table 1 were considered. The corresponding estimates of A(x,y) and ϕ(x,y) reconstructed by the optimal estimation method are shown in Fig. 12. Subfigures (a)–(c) contain the optimal estimates of A(x,y) for detector position error levels (1), (2), and (3), respectively. The corresponding estimates of ϕ(x,y) are contained in subfigures (d)–(f). Despite the geometry uncertainties, the reconstructed A(x,y) and ϕ(x,y) closely resemble the corresponding images in Fig. 8(d) and (h), which were reconstructed from the same noisy data but assuming perfect knowledge of the imaging geometry. The mean square error (MSE) values of the estimates of A(x,y) displayed in Fig. 12(a)–(c) are 1.26×10^-7, 1.27×10^-7, and 1.32×10^-7, respectively. The MSE values for ϕ(x,y) shown in Fig. 12(d)–(f) are 3.37×10^-4, 9.34×10^-4, and 3.37×10^-4, respectively. As a reference, the MSE values for the estimates of A(x,y) and ϕ(x,y) shown in Figs. 8(d) and (h) were 1.25×10^-7 and 4.65×10^-4, respectively. Note that these MSE values were computed by assuming a discrete (Cartesian grid) representation of the true functions A(x,y) and ϕ(x,y), and do not reflect an average over the distribution of the measurement noise. These observations suggest that our multi-plane reconstruction method may be robust to geometry errors under certain conditions. However, a detailed investigation of the effects of geometry errors on the proposed image reconstruction methods remains an important topic for future research.

 

Fig. 8. Estimates of A(x, y) reconstructed from noisy intensity data measured in Geometry ‘A’ by use of detector-planes (a)(1,2), (b)(1,3), (c)(2,3), and (d) an optimally-weighted combination of all three detector-planes. The corresponding estimates of ϕ(x, y) are shown in subfigures (e)–(h).

Download Full Size | PDF

Table 1. Error levels in the detector positions in Geometry ‘A’.

View Table

 

Fig. 9. Empirical variance profiles of Ã(u, v) and ϕ̃(u, v) measured in Geometry ‘A’ are displayed in subfigures (a)–(b), respectively. Each subfigure contains profiles corresponding to the Fourier variances estimated by use of detector-planes (1,2) (dashed curve), detector-planes (1,3) (dashdotted curve), detector-planes (2,3) (dotted curve), and the optimal one (solid curve). Subfigures (c) and (d) display empirical variance profiles of the corresponding images in the spatial domain.

Download Full Size | PDF

7. Summary and conclusions

From knowledge of independent intensity measurements (i.e., phase-contrast radiographs) acquired at distinct states of the imaging system, quantitative phase-contrast X-ray imaging methods seek to reconstruct the projected X-ray absorption and refractive index distributions of an object. If the interaction of the X-ray wavefield with the imaging system is described as a linear, shift-invariant, coherent imaging system, Fourier-based reconstruction formulas can be derived from knowledge of the imaging system’s optical transfer function. A common feature of these reconstruction formulas, however, is the presence of isolated Fourier domain singularities. These poles can greatly amplify noise levels in the estimated Fourier components, which can lead to noisy and/or distorted images in the spatial domain. If intensity measurements are acquired at three or more distinct states of the imaging system, it is sometimes possible to mitigate the noise amplification due to the poles. However, there remains a need for the development of statistically principled reconstruction methods to achieve this.

In this work, we proposed and investigated amethodology for image reconstruction in quantitative phase-contrast imaging when multiple (>2) intensity measurements are available. Linear estimators were proposed that combine the available intensity measurements in such a way that the Fourier components of the desired object properties have reduced variances. These estimators involve the use of combination coefficients that should be chosen to cancel poles present in estimates obtained by use of any two measurements. From knowledge of the measurement noise model, optimal forms of the combination coefficients are derived that minimize the Fourier variance of the final estimates. When such information is not available, we demonstrated that the large noise amplification due to the poles can still be mitigated by an appropriate heuristic choice of the combination coefficients.

 

Fig. 10. Estimates of A(x, y) reconstructed from noisy intensity data measured in Geometry ‘B’ by use of detector-planes (a)(1,2), (b)(1,3), (c)(2,3), and (d) an optimally-weighted combination of all three detector-planes. The corresponding estimates of ϕ(x, y) are shown in subfigures (e)–(h).

Download Full Size | PDF

 

Fig. 11. Empirical variance profiles of Ã(u, v) and ϕ̃(u, v) measured in Geometry ‘B’ are displayed in subfigures (a)–(b), respectively. Each subfigure contains profiles corresponding to the Fourier variances estimated by use of detector-planes (1,2) (dashed curve), detectorplanes (1,3) (dashdotted curve), detector-planes (2,3) (dotted curve), and the optimal one (solid curve). Subfigures (c) and (d) display empirical variance profiles of the corresponding images in the spatial domain.

Download Full Size | PDF

 

Fig. 12. The optimally determined estimates of A(x, y) are reconstructed from noisy intensity data measured in Geometry ‘A’ with detector position uncertainty of (a) error level 1, (b) error level 2, and (c) error level 3. The corresponding estimates of ϕ(x, y) are contained in subfigures (d)–(f).

Download Full Size | PDF

The developed reconstruction methods were investigated, in detail, within the context of propagation-based X-ray phase-contrast imaging. Explicit forms of the optimal estimators were derived in both their continuous and discrete forms. Preliminary computer-simulation studies were conducted to demonstrate the efficacy of the proposed reconstruction methods, and corroborate our theoretical analysis.

As quantitative X-ray phase-contrast imaging is in its infancy, there remain numerous important topics for future investigation. One important topic is the experimental and theoretical investigation of the proposed reconstruction methods within the context of clinically feasible X-ray sources for phase-contrast imaging. This will require generalization of the methods to compensate for non-ideal physical factors such as partial coherence effects. The maximum noise level for which a useful image can be reconstructed by use of the proposed methods is not easily answered by a single number or rule. It will depend, in a non-trivial way, on the measurement geometry, which determines the locations of poles in the reconstruction formulas, the explicit nature of the object’s refractive index distribution, and the ultimate use of the image. It will be important to conduct a detailed investigation of the statistical properties of the reconstructed images and their influence on various task-based detectability measures.