OSA's Digital Library

Virtual Journal for Biomedical Optics

Virtual Journal for Biomedical Optics

| EXPLORING THE INTERFACE OF LIGHT AND BIOMEDICINE

  • Editors: Andrew Dunn and Anthony Durkin
  • Vol. 7, Iss. 1 — Jan. 4, 2012
« Show journal navigation

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

Cheng-Ying Chou and Pin-Yu Huang  »View Author Affiliations


Optics Express, Vol. 19, Issue 24, pp. 24396-24410 (2011)
http://dx.doi.org/10.1364/OE.19.024396


View Full Text Article

Acrobat PDF (1061 KB)





Browse Journals / Lookup Meetings

Browse by Journal and Year


   


Lookup Conference Papers

Close Browse Journals / Lookup Meetings

Article Tools

Share
Citations

Abstract

X-ray phase-contrast tomography (PCT) methods seek to quantitatively reconstruct separate images that depict an object’s absorption and refractive contrasts. Most PCT reconstruction algorithms generally operate by explicitly or implicitly performing the decoupling of the projected absorption and phase properties at each tomographic view angle by use of a phase-retrieval formula. However, the presence of zero-frequency singularity in the Fourier-based phase retrieval formulas will lead to a strong noise amplification in the projection estimate and the subsequent refractive image obtained using conventional algorithms like filtered backprojection (FBP). Tomographic reconstruction by use of statistical methods can account for the noise model and a priori information, and thereby can produce images with better quality over conventional filtered backprojection algorithms. In this work, we demonstrate an iterative image reconstruction method that exploits the second-order statistical properties of the projection data can mitigate noise amplification in PCT. The autocovariance function of the reconstructed refractive images was empirically computed and shows smaller and shorter noise correlation compared to those obtained using the FBP and unweighted penalized least-squares methods. Concepts from statistical decision theory are applied to demonstrate that the statistical properties of images produced by our method can improve signal detectability.

© 2011 OSA

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

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).

3

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

]. 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

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]

6

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).

], and has been exploited to produce three-dimensional (3D) distribution [1

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).

, 7

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

9

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]

] 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

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]

12

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]

]. 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

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]

15

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]

]. 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

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]

, 16

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

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]

]. 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

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

From knowledge of the measured data, the data function Km(x, yr;θ) can be obtained readily as
Km(x,yr;θ)=Im(x,yr;θ)Ii(x,yr;θ)1,
(6)
where Ii(x, yr;θ) is the intensity of the incident wavefield.

3. Image reconstruction

3.1. Phase retrieval

Let (u,vr;θ) denote Ãm,n(u, vr;θ) or ϕ̃m,n(u, vr;θ). If Ãm,n(u, vr;θ) and ϕ̃m,n(u, vr;θ) are assumed to be quasi-bandlimited, containing a maximum frequency Wmax=12Δd, the projection data can be approximated by applying inverse Fourier transform to Ã(u,vr;θ) or ϕ̃(u,vr;θ) described by Eq. (11) or (12) as
P[r,sr;k]1L2p=N2N21qr=N2N21P˜(u=pL,vr=qrL;θ=kΔθ)exp[j2πN(pr+qrsr)].
(13)

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
fa=n=0N21αnψn,
(14)
where αn is the n-th component of the coefficient vector α, {ψn(r)}n=0N21 is a set of expansion functions, and N2 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=Vd3rψn(r)f(r),
(15)
where f(r⃗) describes the 3D refractive index distribution. Although the conventional voxel expansion function is adopted here, other sets of expansion functions [24

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

, 25

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]

] 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α,
(16)
where g ∈ ℝM corresponds to the lexicographically ordered collection of projection estimates, α ∈ ℝN2 denotes the object function, and R ∈ ℝM×N2 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, yr = srΔd) and ϕ(x = rΔd, yr = srΔ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

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]

]. 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

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]

],
Φ(α)=12(gRα)W(gRα)+ηU(α),
(17)
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(α)=12n=0N21k𝒩n(αnαk)2,
(18)
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

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]

]. This will yield an estimate of α̂ that satisfies
α^=argminαΦ(α).
(19)
In phase-contrast tomography, the elements of g refer to the finite-dimensional representation of estimated λ/2π · Am,n[r, sr] or −λ/2π · ϕm,n[r, sr], and α of those refers to the corresponding βm,n[r, s,t] or am,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 Im[r, sr] 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

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]

]
Im[r,sr]=Im0[r,sr]+nm[r,sr],
(20)
where Im0[r,sr] and nm[r, sr] denote the noiseless intensity data and an additive noise term, respectively. We assume that the noise satisfies
Cov{nm[r,sr],nm[r,sr]}=Var{nm[r,sr]}δrrδssrδmm,
(21)
where δ denotes the Kronecker delta function, and Cov{nm[r, sr], nm [r′,sr]} and Var{nm[r, sr]} denote the autocovariance and variance functions of the noise, respectively [28

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

, 29

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).

].

The autocovariance functions of the projected phase and absorption estimates were computed analytically according to Eqs. (31) and (33) in Ref. [14

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]

], which take on the following forms
Cov{ϕm,n(x,yr),ϕm,n(x,yr)}|x=rΔd,yr=srΔdx=rΔd,yr=srΔdCov{ϕm,n[r,sr],ϕm,n[r,sr]}=1N4p=0N1qr=0N1exp[j2πN(pr+qrsr)]×p=0N1qr=0N1exp[j2πN(pr+qrsr)]Cov{ϕ˜m,n[p,qr],ϕ˜m,n[p,qr]},
(22)
and
Cov{Am,n(x,yr),Am,n(x,yr)}|x=rΔd,yr=srΔdx=rΔd,yr=srΔdCov{Am,n[r,sr],Am,n[r,sr]}=1N4p=0N1qr=0N1exp[j2πN(pr+qrsr)]×p=0N1qr=0N1exp[j2πN(pr+qrsr)]Cov{A˜m,n[p,qr],A˜m,n[p,qr]}.
(23)

4.2. Weighting matrix

The weighting function in Eq. (17) was computed by inverting the autocovariance matrix of g, Kg. Because the additive noise is assumed to be uncorrelated between different measurement planes, the autocovariance Kg (M × M) for a transverse slice of the projection estimate is a block diagonal matrix
Kg=(K1000K2000KNv)
(24)
where M = N × Nv and each block size is N × N, with N and Nv denoting the linear dimension of a square detector and the number of tomographic view angles, respectively. Here K1KNv correspond to λ2/4π2 · Cov {ϕm,n(x, yr),ϕm,n(x′, yr)} or λ2/4π2 · Cov{Am,n(x, yr), Am,n(x′, yr)} at the 1-st to Nv-th tomographic projections.

Although the dimension of Kg 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:
Kg1=(K1000K2000KNv)1=(K11000K21000KNv1).
(25)

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 (HWHW) of Hessian matrix that exploits the second-order statistical properties of the projection data corresponding to
HW=RWR,
(26a)
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=RR.
(26b)
In Fig. 2, the normalized eigen spectra of the Hermitian operators HWHW and HH (η = 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 HW [29

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

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

]. 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.

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 RWR and RR, respectively.

The wavelength of the monochromatic incident X-ray beam was 0.8265 Å. The intensity data were acquired on two distinct detector planes positioned at z1 = 9 and z2 = 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, vr) and ϕ̃1,2(u, vr) were reconstructed by use of Eqs. (11) and (12), respectively. From these Fourier data, estimates of A1,2(x, yr) and ϕ1,2(x, yr) 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 π).

5.3. Reconstructed results

Because the presence of singularity in the reconstruction formula (Eq. (12)), it was shown that ϕ(x,yr) is contaminated by low frequency noise while A(x, yr) is not [14

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]

]. 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

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]

]. 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 RWR and RR. 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 RWR and RR, respectively.

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

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]

]. 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 I1[r, sr, k] and I2[r, sr, 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{α^n,α^n}=1J1[i=1Jα^niα^ni1Ji=1Jα^nii=1Jα^ni],
(27)
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{a1,2[0, s,t],a1,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 RWR and RR, respectively.
Fig. 6 Autocovariance profiles Cov{a1,2[0, s,0], a1,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 RWR and RR, respectively.

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{a1,2[0, s,0], a1,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 RWR and RR, respectively.

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α^1[α¯1α¯0], and 1 and 0 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.

The ROC curves that describe the signal detection performance involving the images a1,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 RWR and RR, respectively.

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]

OCIS Codes
(100.5070) Image processing : Phase retrieval
(110.4280) Imaging systems : Noise in imaging systems
(110.7440) Imaging systems : X-ray imaging
(170.3010) Medical optics and biotechnology : Image reconstruction techniques

ToC Category:
Imaging Systems

History
Original Manuscript: September 7, 2011
Revised Manuscript: October 15, 2011
Manuscript Accepted: October 23, 2011
Published: November 14, 2011

Virtual Issues
Vol. 7, Iss. 1 Virtual Journal for Biomedical Optics

Citation
Cheng-Ying Chou and Pin-Yu Huang, "Image reconstruction in phase-contrast tomography exploiting the second-order statistical properties of the projection data," Opt. Express 19, 24396-24410 (2011)
http://www.opticsinfobase.org/vjbo/abstract.cfm?URI=oe-19-24-24396


Sort:  Author  |  Year  |  Journal  |  Reset  

References

  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. SPIE4503, 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. Express11, 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. A19, 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. A23, 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. Express17, 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,” Optik49, 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. SPIE4503, 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]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.


« Previous Article  |  Next Article »

OSA is a member of CrossRef.

CrossCheck Deposited