Image reconstruction in phase-contrast tomography exploiting the second-order statistical properties of the projection data

Cheng-Ying Chou; Pin-Yu Huang

doi:10.1364/OE.19.024396

1. Introduction

X-ray phase-contrast tomography (PCT) is an imaging technique that can produce two separate images that respectively describe an object’s attenuation and phase shift characteristics with respect to the transmitted X-ray wavefield [1–3]. Different from conventional X-ray computed tomography (CT), PCT can exploit variations in the X-ray refractive index distribution of an object. Because variations in the real valued refractive index are several orders of magnitude larger than the imaginary component, it can permit for the visualization of objects with similar or identical absorption contrast. Additionally, phase contrast can persist at higher X-ray energies while absorption contrast decays rapidly, and thus PCT has the potential to reduce radiation dose. Its other attractive features include better stability with respect to high frequency noise and improved spatial locality. Hence it is particularly suitable for biomedical imaging applications [4–6], and has been exploited to produce three-dimensional (3D) distribution [1, 7–9] of an object in phase-contrast tomography.

Quantitative PCT seeks to reconstruct two separate 3D images describing the real and imaginary components of the complex refractive index of objects. In this article, we consider the propagation-based (i.e., in-line) PCT variant. The implementation of the propagation-based PCT is straightforward because no specialized optical elements are needed. Additionally, its tolerance to beam polychromaticity enables the implementation on a benchtop system. This is accomplished by moving the detector away from the contact plane and recording the tomographic scans at different object-to-detector spacings. The tomographic reconstruction process can generally be divided into two steps. First, a phase retrieval formula is applied to yield a collection of projected phase and absorption estimates [10–12]. Subsequently, the 3D distribution of the refractive index can be reconstructed on a slice by slice basis by inverting two-dimensional (2D) Radon transform.

It is noteworthy that the Fourier-based phase retrieval formulas based on transport of intensity equation (TIE) or the contrast transfer function (CTF) contain Fourier domain singularities. The singularities at zero frequency will result in an amplification of low frequency noise in reconstructed phase images [12–15]. The application of the ramp filter in the filtered backprojection (FBP) algorithm can slightly lessen such noisy appearance in the resulting tomographic images. However, the presence of lumpy background can still obscure the detection of signals and affect the accuracy of diagnosis. If the reconstruction algorithm can fully take into account the statistical properties of sinograms, the low frequency noise can be mitigated and the quality of the reconstructed image can be enhanced. Recently, the analytic expressions for image covariance in planar-mode phase-contrast imaging have been determined [14, 16, 17]. The planar-mode images represent the input projection data for PCT. The statistical properties of the projected estimates can be accounted for in the tomographic reconstruction to produce images with better quality over the conventional FBP method and the unweighted penalized least-squares algorithm, which will be abbreviated as PLS hereafter. In this work, we employ a penalized weighted least-squares (PWLS) method that takes into account the autocovariance function of the projected phase images, and demonstrate that the second-order statistical properties of the reconstructed images are favorable and can improve signal detectability.

This paper is organized as follows. The imaging physics of propagation-based PCT is briefly reviewed in Sec. 2. In Sec. 3, the reconstruction formulas employed in this work are introduced. The noise model assumed in the measurements and the weighting function used in the PWLS algorithm are described in Sec. 4. A computer-simulation study was carried out to investigate and evaluate the reconstruction method, and the numerical results of which are presented in Sec. 5. Finally we conclude the article with a discussion and summary in Sec. 6.

2. Imaging physics of propagation-based X-ray phase-contrast tomography

In this section, we briefly review the imaging model of the propagation-based phase-contrast tomography developed by Cloetens [18]. A schematic of the scanning geometry for PCT is illustrated in Fig. 1. To describe the tomographic scanning geometry, we employ the rotated coordinates (x,y_r, z_r) that are related to the reference system (x, y,z) as y_r = y cosθ + zsinθ and z_r = −ysinθ + zcosθ. Here θ denotes the tomographic view angle measured from positive y-axis. Consider that a monochromatic plane wave U_i(x, y_r;θ) with wavelength λ propagates along the z_r-axis and irradiates on the object centered at the origin of the reference coordinate system. The object is characterized by its complex X-ray refractive index distribution as

n (\vec{r}) = 1 - a (\vec{r}) + j β (\vec{r}),

where r⃗ ≡ (x,y,z),

j \equiv \sqrt{- 1}

, and a(r⃗) and β(r⃗) correspond to refractive and absorption properties of the object, respectively. The transmitted wavefield U_t (x, y_r, z_r = 0) on the contact plane immediately behind the object can be described by

U_{t} (x, y_{r}; z_{r} = 0, θ) = U_{i} (x, y_{r}; θ) exp [- A (x, y_{r}; θ) + j ϕ (x, y_{r}; θ)],

where the absorption and phase perturbations are respectively related to the complex refractive index as

A (x, y_{r}) = (2 π / λ) \int_{L} β (\vec{r}) d z_{r}

and

ϕ (x, y_{r}) = - (2 π / λ) \int_{L} a (\vec{r}) d z_{r},

with L denoting the path of X-ray. Equations (3) and (4) are essentially stacks of 2D Radon transforms along x-axis. The intensity of the transmitted wavefield is recorded on two or more detector planes specified by the rotated coordinates (x,y_r) positioned at different z_r > 0. On the detector plane z_r = z_m behind the object, the measured intensity is related to the transmitted wavefield as

I_{m} (x, y_{r}; θ) = {| U_{t} (x, y_{r}; z_{r} = 0, θ) * * h_{z_{m}} (x, y_{r}) |}^{2},

where h_{z_m} (x,y_r) is Fresnel propagator and ‘**’ denotes the 2D convolution operation. The tomographic scanning is performed by keeping the object fixed while simultaneously rotating the X-ray source and the detector, or equivalently, by keeping the imaging system fixed and rotating the object.

Fig. 1 A schematic of the imaging geometry of propagation-based X-ray phase-contrast tomography.

Download Full Size | PDF

From knowledge of the measured data, the data function K_m(x, y_r;θ) can be obtained readily as

K_{m} (x, y_{r}; θ) = \frac{I_{m} (x, y_{r}; θ)}{I_{i} (x, y_{r}; θ)} - 1,

where I_i(x, y_r;θ) is the intensity of the incident wavefield.

3. Image reconstruction

3.1. Phase retrieval

In practical implementation, a digital detector system is often utilized, and as such the measured intensity data will correspond to an ordered collection of numbers instead of a function of a continuous variable. The discrete version of intensity data on a measurement plane z_r = z_m is denoted as

{I_{m} [r, s_{r}; k] = I (x, y_{r}; z_{r} = z_{m}, θ) |}_{x = r Δ_{d}, y_{r} = s_{r} Δ_{d}, θ = k Δ_{θ}},

where r and s_r are integer-valued detector element indices and k is the tomographic view index. Here,

Δ_{d} = \frac{L}{N}

denotes the detector element dimension in a square detector array of dimension L × L and N denotes the number of samples measured in each dimension. Without loss of generality, we will assume that N is an even number. The quantity Δ_θ denotes the angular sampling interval between the uniformly distributed view angles. Equation (7) assumes an idealized (Dirac delta) sampling; namely, the averaging effects of sampling aperture are not considered. However, the analysis that follows can be generalized to address such effects. Hereafter, a square bracket, ‘[·]’, is employed to describe discretely sampled quantities, in order to distinguish them from functions of continuous variables. Let

{\tilde{K}}_{m} [p, q_{r}; k] = \sum_{r = 0}^{N - 1} \sum_{s_{r} = 0}^{N - 1} K_{m} [r, s_{r}; k] exp [- j \frac{2 π}{N} (p r + q_{r} s_{r})]

denote the 2D discrete Fourier transform (DFT) of the discrete data function

K_{m} [r, s_{r}; k] = I_{m} [r, s_{r}; k] / I_{i} - 1,

where [p,q_r] are integer-valued spatial frequency indices conjugated to the detector element indices [r,s_r], and we assume the incident intensity is homogeneous across the detector plane. When the detector sampling interval Δ_d is sufficiently small, so that the effects of aliasing are not large, the continuous and discrete FTs of the data function are related as [19]

{{\tilde{K}}_{m} (u, v_{r}; θ) |}_{u = \frac{p}{L}, v_{r} = \frac{q r}{L}, θ = k Δ_{θ}} \approx Δ_{d}^{2} {\tilde{K}}_{m} [p, q_{r}; k] .

If the imaged object is weakly absorbing (i.e., |A(x,y_r)| << 1) and the variation in ϕ(x,y_r) is smooth, it has been demonstrated that A(x,y_r) and ϕ(x,y_r) can be determined readily by the algebraic formulas in Fourier domain [20–22]. Based on the CTF formulation of phase-contrast imaging [12, 18], Ã(u,v_r;θ) and ϕ̃(u,v_r;θ) can be determined from the data functions as

\begin{array}{l} {\tilde{A}}_{m, n} (u = \frac{p}{L}, v_{r} = \frac{q_{r}}{L}; θ = k Δ_{θ}) \\ = \frac{Δ_{d}^{2}}{D_{m, n} (\frac{p}{L}, \frac{q_{r}}{L})} [sin (π λ z_{n} f_{d}^{2}) {\tilde{K}}_{m} [p, q_{r}; k] - sin (π λ z_{m} f_{d}^{2}) {\tilde{K}}_{n} [p, q_{r}; k]] \end{array}

and

\begin{array}{l} {\tilde{ϕ}}_{m, n} (u = \frac{p}{L}, v_{r} = \frac{q_{r}}{L}; θ = k Δ_{θ}) \\ = \frac{Δ_{d}^{2}}{D_{m, n} (\frac{p}{L}, \frac{q_{r}}{L})} [cos (π λ z_{n} f_{d}^{2}) {\tilde{K}}_{m} [p, q_{r}; k] - cos (π λ z_{m} f_{d}^{2}) {\tilde{K}}_{n} [p, q_{r}; k]], \end{array}

respectively [18, 23], where

f_{d}^{2} \equiv (p^{2} + q_{r}^{2}) \frac{1}{L^{2}}

,

D_{m, n} (\frac{p}{L}, \frac{q_{r}}{L}) \equiv 2 sin (π λ f_{d}^{2} (z_{m} - z_{n}))

, and subscripts (m,n) have been added to Ã(u,v_r;θ) and ϕ̃(u,v_r;θ) to denote that they are estimated by use of measurements K_m[r, s_r; k] and K_n[r, s_r; k]. Note that a singularity is present at the zero frequency u = v_r = 0 in Eq. (12), indicating that the low-frequency components of ϕ will contain greatly amplified noise levels. The pole will affect the noise properties of the reconstructed phase images ϕ and a.

Let P̃(u,v_r;θ) denote Ã_m,n(u, v_r;θ) or ϕ̃_m,n(u, v_r;θ). If Ã_m,n(u, v_r;θ) and ϕ̃_m,n(u, v_r;θ) are assumed to be quasi-bandlimited, containing a maximum frequency $W_{max} = \frac{1}{2 Δ_{d}}$ , the projection data can be approximated by applying inverse Fourier transform to Ã(u,v_r;θ) or ϕ̃(u,v_r;θ) described by Eq. (11) or (12) as

P [r, s_{r}; k] \approx \frac{1}{L^{2}} \sum_{p = - \frac{N}{2}}^{\frac{N}{2} - 1} \sum_{q_{r} = - \frac{N}{2}}^{\frac{N}{2} - 1} \tilde{P} (u = \frac{p}{L}, v_{r} = \frac{q_{r}}{L}; θ = k Δ_{θ}) exp [j \frac{2 π}{N} (p r + q_{r} s_{r})] .

3.2. Tomographic reconstruction

Following the computation of projection data, the tomographic reconstruction was carried out. Because we focus on iterative reconstruction algorithms here, we will also require a discrete representation of the object

f_{a} = \sum_{n = 0}^{N^{2} - 1} α_{n} ψ_{n},

where α_n is the n-th component of the coefficient vector α,

{ψ_{n} (\vec{r})}_{n = 0}^{N^{2} - 1}

is a set of expansion functions, and N² corresponds to the number of expansion functions. In this work, we adopt the conventional voxel expansion as the choice for ψ_n, and the corresponding coefficient component is given by

α_{n} = \int_{V} d^{3} \vec{r} ψ_{n} (\vec{r}) f (\vec{r}),

where f(r⃗) describes the 3D refractive index distribution. Although the conventional voxel expansion function is adopted here, other sets of expansion functions [24, 25] could be employed to form a finite-dimensional approximate object representation.

The tomographic reconstruction problem can be approximated by the discrete linear model

g = R α,

where g ∈ ℝ^M corresponds to the lexicographically ordered collection of projection estimates, α ∈ ℝ^N² denotes the object function, and R ∈ ℝ^M×N² represents the discrete approximation of 2D Radon transform. The application of Eqs. (11) and (12) followed by the inverse DFT in Eq. (13) will yield the estimation of A(x = rΔ_d, y_r = s_rΔ_d) and ϕ(x = rΔ_d, y_r = s_rΔ_d) at each view angle, respectively. Subsequently, β[r,s,t] and a[r,s,t] can be obtained by inverting 2D Radon transforms. This can be accomplished by applying a statistical reconstruction method aimed at finding the solution that minimizes the mismatch between the observed data and that projected by the estimated image. Unregularized reconstruction methods produce images of increasingly noisy appearance with iteration [26]. Incorporating smoothness penalty or prior information can improve the stability of the tomographic reconstruction.

A PWLS method that takes into consideration of the statistical properties of the observed data is an useful approach to obtain the solution. The cost function Φ(α) that we aim to minimize takes on the following form [27],

Φ (α) = \frac{1}{2} {(g - R α)}^{†} W (g - R α) + η U (α),

where W corresponds to the inverse of the autocovariance matrix of g, † denotes adjoint operation, U(α) is the penalty function that smoothes the object function α, and η is the regularization parameter that controls the tradeoff between the noise and resolution. In this study we employed a simple quadratic smoothness penalty given by

U (α) = \frac{1}{2} \sum_{n = 0}^{N^{2} - 1} \sum_{k \in 𝒩_{n}} {(α_{n} - α_{k})}^{2},

where 𝒩_n is the set of 4 neighbors of the n-th pixel. The objective of the reconstruction method is to find an image that minimizes the cost function. This can be accomplished by employing optimization algorithms such as conjugate-gradient (CG) methods [27]. This will yield an estimate of α̂ that satisfies

\hat{α} = arg min_{α} Φ (α) .

In phase-contrast tomography, the elements of g refer to the finite-dimensional representation of estimated λ/2π · A_m,n[r, s_r] or −λ/2π · ϕ_m,n[r, s_r], and α of those refers to the corresponding β_m,n[r, s,t] or a_m,n[r, s,t].

4. Noise model

4.1. Noise statistics of measured intensity and reconstructed images

In the presence of stochastic noise, the measured intensity I_m[r, s_r] represents a stochastic quantity that contains the noiseless intensity plus noise. In order not to further complicate the noise analysis of the proposed iterative method, without loss of generality, we considered a simplified measurement model [13]

I_{m} [r, s_{r}] = I_{m}^{0} [r, s_{r}] + n_{m} [r, s_{r}],

where

I_{m}^{0} [r, s_{r}]

and n_m[r, s_r] denote the noiseless intensity data and an additive noise term, respectively. We assume that the noise satisfies

Cov {n_{m} [r, s_{r}], n_{m^{'}} [r^{'}, s_{r}^{'}]} = Var {n_{m} [r, s_{r}]} δ_{r r^{'}} δ_{s s_{r}^{'}} δ_{m m^{'}},

where δ denotes the Kronecker delta function, and Cov{n_m[r, s_r], n_m′ [r′,s′_r]} and Var{n_m[r, s_r]} denote the autocovariance and variance functions of the noise, respectively [28, 29].

The autocovariance functions of the projected phase and absorption estimates were computed analytically according to Eqs. (31) and (33) in Ref. [14], which take on the following forms

\begin{array}{l} {Cov {ϕ_{m, n} (x, y_{r}), ϕ_{m, n} (x^{'}, y_{r}^{'})} |}_{\begin{array}{l} x = r Δ_{d}, y_{r} = s_{r} Δ_{d} \\ x^{'} = r^{'} Δ_{d}, y_{r}^{'} = s_{r}^{'} Δ_{d} \end{array}} \approx Cov {ϕ_{m, n} [r, s_{r}], ϕ_{m, n} [r^{'}, s_{r}^{'}]} \\ = \frac{1}{N^{4}} \sum_{p = 0}^{N - 1} \sum_{q_{r} = 0}^{N - 1} exp [j \frac{2 π}{N} (p r + q_{r} s_{r})] \\ \times \sum_{p^{'} = 0}^{N - 1} \sum_{q_{r}^{'} = 0}^{N - 1} exp [- j \frac{2 π}{N} (p^{'} r^{'} + q_{r}^{'} s_{r}^{'})] Cov {{\tilde{ϕ}}_{m, n} [p, q_{r}], {\tilde{ϕ}}_{m, n} [p^{'}, q_{r}^{'}]}, \end{array}

and

\begin{array}{l} {Cov {A_{m, n} (x, y_{r}), A_{m, n} (x^{'}, y_{r}^{'})} |}_{\begin{array}{l} x = r Δ_{d}, y_{r} = s_{r} Δ_{d} \\ x^{'} = r^{'} Δ_{d}, y_{r}^{'} = s_{r}^{'} Δ_{d} \end{array}} \approx Cov {A_{m, n} [r, s_{r}], A_{m, n} [r^{'}, s_{r}^{'}]} \\ = \frac{1}{N^{4}} \sum_{p = 0}^{N - 1} \sum_{q_{r} = 0}^{N - 1} exp [j \frac{2 π}{N} (p r + q_{r} s_{r})] \\ \times \sum_{p^{'} = 0}^{N - 1} \sum_{q_{r}^{'} = 0}^{N - 1} exp [- j \frac{2 π}{N} (p^{'} r^{'} + q_{r}^{'} s_{r}^{'})] Cov {{\tilde{A}}_{m, n} [p, q_{r}], {\tilde{A}}_{m, n} [p^{'}, q_{r}^{'}]} . \end{array}

4.2. Weighting matrix

The weighting function in Eq. (17) was computed by inverting the autocovariance matrix of g, K_g. Because the additive noise is assumed to be uncorrelated between different measurement planes, the autocovariance K_g (M × M) for a transverse slice of the projection estimate is a block diagonal matrix

K_{g} = (\begin{matrix} K_{1} & 0 & \dots & 0 \\ 0 & K_{2} & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & K_{N_{v}} \end{matrix})

where M = N × N_v and each block size is N × N, with N and N_v denoting the linear dimension of a square detector and the number of tomographic view angles, respectively. Here K₁ ⋯ K_{N_v} correspond to λ²/4π² · Cov {ϕ_m,n(x, y_r),ϕ_m,n(x′, y′_r)} or λ²/4π² · Cov{A_m,n(x, y_r), A_m,n(x′, y′_r)} at the 1-st to N_v-th tomographic projections.

Although the dimension of K_g is enormous, its inversion is not intractable. The inverse of a block diagonal matrix is another block diagonal matrix, composed of the inverse of each block, as follows:

K_{g}^{- 1} = {(\begin{matrix} K_{1} & 0 & \dots & 0 \\ 0 & K_{2} & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & K_{N_{v}} \end{matrix})}^{- 1} = (\begin{matrix} K_{1}^{- 1} & 0 & \dots & 0 \\ 0 & K_{2}^{- 1} & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & K_{N_{v}}^{- 1} \end{matrix}) .

The 3D distribution of the refractive index can be reconstructed by applying the proposed PWLS algorithm on a slice by slice basis. Taking example of a 4k × 4k detector, the reconstruction involves the inversion of each block matrix with dimensions of 4000 × 4000, which can be computed readily. Consequently, even though the PWLS method involves computation of matrix inversion, it is applicable to real data of such dimensions.

5. Computer-simulation studies

Computer simulation studies were conducted to investigate how the addition of the autoco-variance matrix in the reconstruction method will influence the statistical properties of reconstructed images in PCT. Subsequent to its empirical determination, the autocovariance function of the tomographic phase image was employed to implement observer studies.

5.1. Eigen analysis of Hessian matrices

In order to analyze the effects of the autocovariance function in the system matrix, we performed the eigen analysis on the Hermitian operator (H_W^†H_W) of Hessian matrix that exploits the second-order statistical properties of the projection data corresponding to

H_{W} = R^{†} W R,

where the weighting matrix W was computed when the standard deviation of Gaussian noise was set as σ = 1%. On the other hand, the system matrix of the unweighted least-squares method was also computed for comparison purpose. The unweighted system matrix is given by

H = R^{†} R .

In Fig. 2, the normalized eigen spectra of the Hermitian operators H_W^†H_W and H^†H (η = 0 for both cases) are respectively depicted by solid and dashed curves. The decay behavior of eigenvalues corresponding to the unweighted Hessian matrix is greatly improved when the inverse autocovariance function is included in H_W [29,30]. Note that the semi-logarithmic axes were employed to show both curves in the same plot. For the noise model considered here, the profiles of normalized eigen spectra remain unchanged with respect to the noise level.

Fig. 2 Normalized eigenspectra of the weighted (solid curve) and unweighted (dashed curve) Hessian matrices.

Download Full Size | PDF

5.2. Numerical phantom and simulated measurement data

We employed two sets of mathematical phantoms to represent the objects of interest. Phantom A comprised of 5 uniform ellipsoids was employed to demonstrate the images reconstructed by the FBP, PWLS and PLS algorithms at low noise level of σ = 1%. The reconstructed images of Phantom A by use of these algorithms are presented in Fig. 3. At higher noise level of σ = 10%, we employed Phantom set B that consists of one small ellipsoid, the signal, embedded in one large ellipsoid, the background, as our numerical phantoms to investigate the second-order statistical properties of the reconstructed images and the associated signal detection studies. Additional description about Phantom set B is provided in Sec. 5.5.

Fig. 3 The images reconstructed by use of the (a) FBP and PWLS algorithms employing (b) 10, (c) 50, (d) 90 iterations. The bottom row contains images reconstructed using the PLS method employing (e) 3, (f) 7, (g) 11, (h) 15 iterations. The standard deviation of the additive Gaussian noise and regularization parameter were set as σ = 1% and η = 0.0434 for the normalized matrices R^†WR and R^†R, respectively.

Download Full Size | PDF

The wavelength of the monochromatic incident X-ray beam was 0.8265 Å. The intensity data were acquired on two distinct detector planes positioned at z₁ = 9 and z₂ = 38 mm behind the object, respectively. The detector was assumed to contain 128 × 128 elements of dimension of 5 μm. Noisy intensity data were produced according to Eq. (20), where the standard deviation of Gaussian noise was σ = 1%. Estimates of the Fourier components of Ã_1,2(u, v_r) and ϕ̃_1,2(u, v_r) were reconstructed by use of Eqs. (11) and (12), respectively. From these Fourier data, estimates of A_1,2(x, y_r) and ϕ_1,2(x, y_r) were obtained after application of the 2D inverse Fourier transform. A set of tomographic data was obtained by repeating the above procedures at 180 evenly distributed tomographic view angles θ over [0 π).

The autocovariance functions of the projected estimates were computed analytically by use of Eqs. (22) and (23) (i.e., Eqs. (31) and (33) in Ref. [14]). The elements of the weighting matrix W in Eq. (17) were specified by inverting the computed autocovariance matrix, according to Eq. (25). In order to search for the solution that minimizes the objective function expressed in Eq. (17), the conjugate-gradient method was employed. In this work, we used the Polak-Ribiere CG method to calculate the search direction and solve the inverse problem iteratively [27, 31].

5.3. Reconstructed results

Because the presence of singularity in the reconstruction formula (Eq. (12)), it was shown that ϕ(x,y_r) is contaminated by low frequency noise while A(x, y_r) is not [14]. Next, the collection of phase estimates was used as the projection data for tomographic reconstruction. This was accomplished by use of FBP, PWLS and PLS algorithms. Although the ramp filter in the classical FBP algorithm can help to mitigate the Fourier pole at zero frequency, the reconstructed refractive index images still possess strong low-frequency noise, as compared to the attenuation images [15]. By use of the measured intensity data with noise level σ = 1%, the reconstructed estimates of a(r⃗) corresponding to the transverse slice x = 0 are presented in Fig. 3. The regularization parameter was chosen as η = 0.0434 for the normalized matrices R^†WR and R^†R. Subfigures (a) and (b)–(d) of Fig. 3 correspond to the images reconstructed by use of the FBP and PWLS algorithms at the 10-th, 50-th, and 90-th iterations, respectively.

The bottom row contains the images reconstructed by use of the PLS algorithm with sub-figures (e)–(h) corresponding to images produced at the 3-rd, 7-th, 11-th and 15-th iterations, respectively. Note that the reconstructed images exhibit excessive low-frequency noise appearance at early iterations, but later show strong high-frequency noise contamination as iteration increases.

Subsequent to the image reconstruction, the corresponding cost function ||Φ(a^(ℓ))|| and normalized mean square error of reconstructed images versus iteration number defined as ||â^(ℓ) – a||/||a|| were examined and are presented in Fig. 4, where â^(ℓ) is the estimate in the ℓ-th iteration and ||·|| denotes L2 norm. The normalized mean square errors in Figs. 4(b) and 4(d) reveal that the PWLS and PLS algorithms converge differently. The proposed PWLS method, though converges slower than PLS one, reaches the solution of minimal error around the 80-th iteration. As for the PLS algorithm, whose minimal error solution occurs around the 3-rd iteration diverges soon after that. Based on the above observation, in all the simulation studies hereafter, we terminated the CG algorithm after 90 and 15 iterations for the PWLS and PLS algorithms, respectively.

Fig. 4 The cost function and normalized mean square error at each iteration for the PWLS method are contained in subfigures (a) and (b). The corresponding cost function and normalized mean square error for the PLS method are contained in subfigures (c) and (d), respectively. The additive Gaussian noise and regularization parameter were set as σ = 1% and η = 0.0434 for the normalized matrices R^†WR and R^†R, respectively.

Download Full Size | PDF

5.4. Empirical determination of reconstructed image statistics

The statistical properties of phase-contrast tomography employing the FBP algorithm were computed analytically using Eq. (35) and the equations in Sec. III C in Ref. [15]. The statistical properties of images produced by our proposed PWLS and the conventional PLS algorithms were computed from an ensemble of J (=2000) reconstructed refractive images estimated up to 90 and 15 iterations from noisy realizations of intensity measurement pairs I₁[r, s_r, k] and I₂[r, s_r, k], as described in Eq. (20) with σ = 1%. J noisy sets of α̂ were reconstructed by the PWLS and PLS methods and their autocovariance functions were empirically determined by

Cov {{\hat{α}}_{n}, {\hat{α}}_{n^{'}}} = \frac{1}{J - 1} [\sum_{i = 1}^{J} {\hat{α}}_{n}^{i} {\hat{α}}_{n^{'}}^{i} - \frac{1}{J} \sum_{i = 1}^{J} {\hat{α}}_{n}^{i} \sum_{i = 1}^{J} {\hat{α}}_{n^{'}}^{i}],

where superscript i denotes the i-th noisy realization.

The autocovariance functions of the refractive images obtained by use of the FBP and PWLS methods are contained in Fig. 5(a), and Figs. 5(b)–5(d), respectively. The autocovariance functions corresponding to those images in Figs. 3(e)–(h) are shown in Figs. 5(e)–(h), indicating that the PLS results are greatly influenced by low frequency noise after a few iterations and also by high frequency noise with the application of more iterations. In contrast to the corresponding FBP result, both the magnitude and length of noise correlation of the refractive images are reduced by use of the PWLS method despite of the noise amplification as iterations proceed. This is confirmed by the autocovariance profiles depicted by Figs. 6(a)–(c), in which PWLS results (solid curves) obtained after (a) 10, (b) 50, and (c) 90 iterations show shorter lengths of noise correlation than the PLS result (solid curve with circle marker) obtained after 3 iterations.

Fig. 5 Top row shows the autocovariance functions Cov{a_1,2[0, s,t],a_1,2[0, 0,0]} reconstructed by use of the (a) FBP and PWLS methods after (b) 10, (c) 50, and (d) 90 iterations. The Gaussian noise and regularization parameter were set as σ = 1% and η = 0.0434 for the normalized matrices R^†WR and R^†R, respectively.

Download Full Size | PDF

Fig. 6 Autocovariance profiles Cov{a_1,2[0, s,0], a_1,2[0, 0,0]} corresponding to refractive images produced by use of the PWLS, PLS and FBP algorithms are denoted by solid, solid with circle marker, and dashed curves, respectively. Solid curves denote the PWLS results estimated after (a) 10, (b) 50, and (c) 90 iterations, while the curve with circle marker for each subfigure denotes the PLS result in the 3-rd iteration. The standard deviation of the Gaussian noise and the regularization parameter were set as σ = 10% and η = 1.4942 for the normalized matrices R^†WR and R^†R, respectively.

Download Full Size | PDF

When the noise level increases to σ = 10%, the noise amplification in the images produced by use of the FBP and PLS methods are much more pronounced than in those produced by use of the PWLS method. The profiles of autocovariance functions of these results are presented in Figs. 7(a) and 7(b), in which the peak values of the profiles for the FBP and PLS results are about 100 times of their counterparts in Fig. 6, while the corresponding noise amplification by use of the PWLS method only increases by about 6-fold.

Fig. 7 Autocovariance profiles Cov{a_1,2[0, s,0], a_1,2[0, 0,0]} corresponding to refractive images generated by use of the (a) FBP (dashed curve) algorithm and PLS (solid curve with circle marker) algorithm at the 3-rd iteration, and by use of the (b) PWLS algorithm after 10 iterations (dash-dotted curve) and after 90 iterations (solid curve), respectively. The standard deviation of the Gaussian noise was σ = 10% and the regularization parameter was η = 1.4942 for the normalized matrices R^†WR and R^†R, respectively.

Download Full Size | PDF

5.5. Signal detection studies

In imaging science, there are many different metrics to evaluate image quality. However, in order to improve diagnostic and/or prognostic accuracy, objective image quality should be assessed based on the intended use of images. In this work, our task of interest is a binary signal detection task (e.g., whether a tumor/lesion is present or not). An observer study can provide a quantitative, reproducible measure of image quality and thus is frequently employed to evaluate technology in use, reconstruction algorithm, and imaging system optimization. But it is too expensive and time consuming to employ a human observer to perform such a repetitive task. Therefore, a model observer can be chosen as an alternative that predicts the performance of a human observer.

Here we employed the Hotelling observer as our model observer, and performed the signal-known-exactly/background-known-exactly (SKE/BKE) detection task by computing the test statistic t(α) = w^†α for 2000 pairs of reconstructed refractive images, where $w = K_{\hat{α}}^{- 1} [{\bar{α}}_{1} - {\bar{α}}_{0}]$ , and ᾱ₁ and ᾱ₀ correspond to the average images of signal present and signal absent, respectively. In the signal detection simulation, Phantom set B was employed and the signal to be detected corresponds to one small ellipsoid that is embedded in another larger ellipsoid (i.e., the background). Here, α refers to the refractive estimates. The numerical phantoms corresponding to signal present and signal absent cases are illustrated in Figs. 8(a) and (b). The figure of merit is the receiver operating characteristic (ROC) curve, which summarizes the performance of the mathematical observer.

Fig. 8 The numerical phantoms corresponding to (a) signal present and (b) signal absent cases.

Download Full Size | PDF

The ROC curves that describe the signal detection performance involving the images a_1,2[0, s,t] reconstructed from noisy measurement data with σ = 10% are contained in Fig. 9, where solid, solid with circle marker and dashed curves denote the results obtained by use of the PWLS, PLS and FBP algorithms, respectively. Note that the solid curve and solid curve with circle marker in Fig. 9 correspond to the PWLS and PLS results obtained at the 50-th and 3-rd iterations. Even though only a single ROC curve corresponding to the PWLS algorithm is shown here, its differences with that evaluated after 70 or 90 iterations are negligible. The ROC curves produced by use of the PWLS algorithm that employed 70 and 90 iterations were found to overlap with the results obtained by use of 50 iterations shown in Fig. 9. The ROC curves indicate that the observer performed considerably better with images reconstructed by use of the PWLS method than those obtained by the conventional PLS and FBP algorithms.

Fig. 9 The ROC curves for refractive images generated by use of the PWLS (solid curve), PLS (solid curve with circle marker) and FBP (dashed curve) algorithms when the standard deviation of the additive Gaussian noise and the regularization parameter were set as σ = 10% and η = 1.4942 for the normalized matrices R^†WR and R^†R, respectively.

Download Full Size | PDF

6. Summary and conclusions

In this work, we employed the PWLS method that utilized the analytically computed autoco-variance function of the projected phase as a priori information in PCT reconstruction. The addition of second-order statistical properties of the projection data not only enhanced the stability of tomographic reconstruction, but also mitigated the low frequency noise in PCT images, as well as the length of noise correlation. This is confirmed by comparing the empirically determined autocovariance profiles of the PWLS results with those estimated by the PLS and FBP algorithms. It is worth noted that the correlation length and magnitude of autocovariance functions increase as iterations proceed, but they both gradually converge when approaching the minimum error solution. The peak value and correlation length of the autocovariance profile of the PWLS results estimated at 90-th iteration are still lower and smaller than those of the FBP and PLS algorithms. The direct impact of such noise reduction in the practical usage of these images is improvement of signal detection, which is corroborated by our numerical observer study.

This study aims to provide some insights into the efficacy of the proposed PWLS method in mitigating low frequency noise in PCT and investigated their effects on the subsequent usage of these images. Notwithstanding the observer studies shown in this work utilized images reconstructed at the 50-th iteration by use of the PWLS algorithm, an earlier or later termination of CG algorithm can be chosen and may result in similar observer performance. In this paper, we considered the parallel incident beam with perfect coherence for the FBP, PLS, and PWLS algorithms. But our analysis can be readily extended to other phase-contrast imaging methods. Although we did not take into account all the physical factors like beam coherence or scattering effects, the study can be readily generalized to address these effects. Other effects like background lumpiness and beam polychromaticity will be the subjects of our future study.

Acknowledgments

This work was supported in part by National Science Council under grant No. NSC 99-2221-E-002-043-MY3.

References and links

1. P. Cloetens, W. Ludwig, E. Boller, L. Helfen, L. Salvo, R. Mache, and M. Schlenker, “Quantitative phase-contrast tomography using coherent synchrotron radiation,” in Developments in X-Ray Tomography III, U. Bonse, ed., Proc. SPIE 4503, 82–91 (2002).

2. S. Mayo, T. Davis, T. Gureyev, P. Miller, D. Paganin, A. Pogany, A. Stevenson, and S. Wilkins, “X-ray phase-contrast microscopy and microtomography,” Opt. Express 11, 2289–2302 (2003). [CrossRef] [PubMed]

3. D. M. Paganin, Coherent X-Ray Optics (Oxford University Press, New York, 2006).

4. R. A. Lewis, “Medical phase contrast x-ray imaging: current status and future prospects,” Phys. Med. Biol. 49, 3573–3583 (2004). URL http://stacks.iop.org/0031-9155/49/3573. [CrossRef] [PubMed]

5. S. Fiedler, A. Bravin, J. Keyrilainen, M. Fernandaz, P. Suortti, W. Thomlinson, M. Tenhenun, P. Virkkunen, and M. Karjalainen-Lindsberg, “Imaging lobular breast carcinoma: comparison of synchrotron radiation CT-DEI technique with clinical CT, mammography and histology,” Phys. Med. Biol. 49, 1–15 (2004). [CrossRef]

6. F. Arfelli, V. Bonvicini, A. Bravin, G. Cantatore, E. Castelli, L. D. Palma, M. D. Michiel, M. Fabrizioli, R. Longo, R. H. Menk, A. Olivo, S. Pani, D. Pontoni, P. Poropat, M. Prest, A. Rashevsky, M. Ratti, L. Rigon, G. Tromba, A. Vacchi, E. Vallazza, and F. Zanconati, “Mammography with synchrotron radiation: phase-detection techniques,” Radiol. 215, 286–293 (2000).

7. A. V. Bronnikov, “Theory of quantitative phase-contrast computed tomography,” J. Opt. Soc. Am. A 19, 472–480 (2002). [CrossRef]

8. T. E. Gureyev, Y. I. Nesterets, D. M. Paganin, and S. W. Wilkins, “Effects of incident illumination on in-line phase-contrast imaging,” J. Opt. Soc. Am. A 23, 34–42 (2006). URL http://josaa.osa.org/abstract.cfm?URI=josaa-23-1-34. [CrossRef]

9. P. McMahon, A. Peele, D. Paterson, K. A. Nugent, A. Snigirev, T. Weitkamp, and C. Rau, “X-ray tomographic imaging of the complex refractive index,” Appl. Phys. Lett. 83, 1480–1482 (2003). [CrossRef]

10. K. A. Nugent, T. E. Gureyev, D. F. Cookson, D. M. Paganin, and Z. Barnea, “Quantitative phase imaging using hard x-rays,” Phys. Rev. Lett. 77, 2961–2964 (1996). [CrossRef] [PubMed]

11. P. Cloetens, M. Pateyron-Salome, J. Y. Buffiere, G. Peix, J. Baruchel, F. Peyrin, and M. Schlenker, “Observation of microstructure and damage in materials by phase sensitive radiography and tomography,” J. Appl. Phys. 81, 5878–5886 (1997). [CrossRef]

12. M. Langer, P. Cloetens, J. P. Guigay, and F. Peyrin, “Quantitative comparison of direct phase retrieval algorithms in in-line phase tomography,” Med. Phys. 35, 4556–66 (2008). [CrossRef] [PubMed]

13. D. Paganin, A. Barty, P. J. McMahon, and K. A. Nugent, “Quantitative phase-amplitude microscopy III. The effects of noise,” J. Microsc. 214, 51–61 (2003). [CrossRef]

14. C.-Y. Chou and M. A. Anastasio, “Influence of imaging geometry on noise texture in quantitative in-line X-ray phase-contrast imaging,” Opt. Express 17, 14,466–14,480 (2009). [CrossRef]

15. C.-Y. Chou and M. A. Anastasio, “Noise texture and signal detectability in propagation-based x-ray phase-contrast tomography,” Med. Phys. 37, 270 (2010). [CrossRef] [PubMed]

16. C.-Y. Chou and M. A. Anastasio, “Influence of imaging geometry on noise texture in x-ray in-line phase-contrast imaging,” in Medical Imaging 2008: Physics of Medical Imaging, J. Hsieh and E. Samei, eds., Proc. SPIE6913, 69131Z (2008). URL http://link.aip.org/link/?PSI/6913/69131Z/1.

17. C.-Y. Chou and M. A. Anastasio, “Statistical properties of X-ray phase-contrast tomography,” in Annual International Conference of the IEEE Engineering in Medicine and Biology Society, 2009. EMBC 2009 , pp. 6648 –6650 (2009). [CrossRef] [PubMed]

18. P. Cloetens, “Contribution to Phase Contrast Imaging, Reconstruction and Tomography with Hard Synchrotron Radiation: Principles, Implementation and Applications,” Ph.D. thesis, Vrije Universiteit Brussel (1999).

19. W. D. Stanley, G. R. Dougherty, and R. Dougherty, Digital Signal Processing, 2nd ed. (Reston Publishing Company, Inc., Reston, VA, 1984).

20. J.-P. Guigay, “Fourier transform analysis of Fresnel diffraction patterns and in-line holograms,” Optik 49, 121–125 (1977).

21. X. Wu and H. Liu, “Clinical implementation of X-ray phase-contrast imaging: theoretical foundations and design considerations,” Med. Phys. 30, 2169–2179 (2003). [CrossRef] [PubMed]

22. C. J. Kotre and I. P. Birch, “Phase contrast enhancement of x-ray mammography: a design study,” Phys. Med. Biol. 44, 2853–2866 (1999). [CrossRef] [PubMed]

23. P. Cloetens, W. Ludwig, E. Boller, L. Helfen, L. Salvo, R. Mache, and M. Schlenker, “Quantitative phase contrast tomography using coherent synchrotron radiation,” in Developments in X-Ray Tomography III, U. Bonse, ed., Proc. SPIE 4503, 82–91 (2001).

24. R. Lewitt, “Alternatives to voxels for image representation in iterative reconstruction algorithms,” Phys. Med. Biol. 37, 705–716 (1992). [CrossRef] [PubMed]

25. M. Defrise and G. T. Gullberg, “Image reconstruction,” Phys. Med. Biol. 51, R139 (2006). URL http://stacks.iop.org/0031-9155/51/i=13/a=R09. [CrossRef] [PubMed]

26. D. Snyder and M. Miller, “The use of sieves to stabilize images produced with the EM algorithm for emission tomography,” IEEE Trans. Nucl. Sci. 32, 3864–3872 (1985). [CrossRef]

27. J. Fessler and S. Booth, “Conjugate-gradient preconditioning methods for shift-variant PET image reconstruction,” IEEE Trans. Image Process. 8, 688–699 (1999). [CrossRef]

28. A. Papoulis and S. U. Pillai, Probability, Random Variables, and Stochastic Processes (McGraw Hill, New York, 2002).

29. H. H. Barrettt and K. J. Myers, Foundations of Image Science, Wiley Series in Pure and Applied Optics (John Wiley & Sons, Inc., Hoboken, New Jersey, 2004).

30. M. Bertero, Introduction to inverse problems in imaging (Taylor & Francis, 1998). [CrossRef]

31. E. Mumcuoglu, R. Leahy, S. R. Cherry, and Z. Zhou, “Fast gradient-based methods for Bayesian reconstruction of transmission and emission PET images,” IEEE Trans. Med. Imag. 13, 687–701 (1994). [CrossRef]

Image reconstruction in phase-contrast tomography exploiting the second-order statistical properties of the projection data

Abstract

1. Introduction

2. Imaging physics of propagation-based X-ray phase-contrast tomography

3. Image reconstruction

3.1. Phase retrieval

3.2. Tomographic reconstruction

4. Noise model

4.1. Noise statistics of measured intensity and reconstructed images

4.2. Weighting matrix

5. Computer-simulation studies

5.1. Eigen analysis of Hessian matrices

5.2. Numerical phantom and simulated measurement data

5.3. Reconstructed results

5.4. Empirical determination of reconstructed image statistics

5.5. Signal detection studies

6. Summary and conclusions

Acknowledgments

References and links

Cited By

Figures (9)

Equations (28)

Optics Express