OSA's Digital Library

Optics Express

Optics Express

  • Editor: C. Martijn de Sterke
  • Vol. 15, Iss. 16 — Aug. 6, 2007
  • pp: 10002–10025
« Show journal navigation

Image reconstruction in quantitative X-ray phase-contrast imaging employing multiple measurements

Cheng-Ying Chou, Yin Huang, Daxin Shi, and Mark A. Anastasio  »View Author Affiliations


Optics Express, Vol. 15, Issue 16, pp. 10002-10025 (2007)
http://dx.doi.org/10.1364/OE.15.010002


View Full Text Article

Acrobat PDF (7916 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 imaging is a technique that aims to reconstruct the projected absorption and refractive index distributions of an object. One common feature of reconstruction formulas for phase-contrast imaging is the presence of isolated Fourier domain singularities, which can greatly amplify the noise levels in the estimated Fourier domain and lead to noisy and/or distorted images in spatial domain. In this article, we develop a statistically optimal reconstruction method that employs multiple (>2) measurement states to mitigate the noise amplification effects due to singularities in the reconstruction formula. Computer-simulation studies are carried out to quantitatively and systematically investigate the developed method, within the context of propagation-based X-ray phase-contrast imaging. The reconstructed images are shown to possess dramatically reduced noise levels and greatly enhanced imaging contrast.

© 2007 Optical Society of America

1. Introduction

Due to the advent and accessibility of synchrotron and laboratory X-ray sources [1

1. A. Krol, A. Ikhlef, J.-C. Kieffer, D. Bassano, C. C. Chamberlain, Z. Jiang, H. Pepin, and S. C. Prasad, “Laser-based microfocused X-ray source for mammography: feasibility study,” Med. Phys. 24, 725–732 (1997). [CrossRef] [PubMed]

, 2

2. R. Waynant, “Toward practical coherent x-ray sources: Potential medical applications,” IEEE J. Quantum Electron. 6, 1465–1469 (2000). [CrossRef]

, 3

3. H. Yamada, “Novel x-ray source based on a tabletop synchrotron and its unique features,” Nucl. Instrum. Methods Phys. Res. B 199, 509–516 (2003).

] that possess high degrees of coherence, new imaging methods, broadly characterized as X-ray phase-contrast imaging methods [4

4. R. Lewis, “Medical phase contrast x-ray imaging: current status and future prospects,” Phys. Med. Biol. 49, 3573–3583 (2004). [CrossRef] [PubMed]

, 5

5. W. Thomlinson, P. Suortti, and D. Chapman, “Recent advances in synchrotron radiation medical research,” Nucl. Instrum. Methods Phys. Res. A 543, 288–296 (2005).

, 6

6. T. Davis, D. Gao, T. E. Gureyev, A. Stevenson, and S. Wilkins, “Phase-contrast imaging of weakly absorbing materials using hard X-rays,” Nature (London) 373, 335–338 (1996).

, 7

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

, 8

8. D. Chapman, W. Thomlinson, R. E. Johnston, D. Washburn, E. Pisano, N. Gmur, Z. Zhong, R. Menk, F. Arfelli, and D. Sayers, “Diffraction enhanced x-ray imaging,” Phys. Med. Biol. 42, 2015–2025 (1997). [CrossRef] [PubMed]

, 9

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

, 10

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

, 11

11. M. N. Wernick, O. Wirjadi, D. Chapman, Z. Zhong, N. P. Galatsanos, Y. Yang, J. G. Brankov, O. Oltulu, M. A. Anastasio, and C. Muehleman, “Multiple-image radiography,” Phys. Med. Biol. 48, 3875–3895 (2003). [CrossRef]

, 12

12. E. F. Donnelly, R. R. Price, and D. R. Pickens, “Characterization of the phase-contrast radiography edge-enhancement effect in a cabinet x-ray system,” Med. Phys. 30, 2292–2296 (2003). [CrossRef] [PubMed]

, 13

13. F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Phase retrieval and differential phase-contrast imaging with low-brilliance X-ray sources,” Nature Phys. 2, 258–261 (2006). [CrossRef]

] 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

14. B.L. Henke, E.M. Gullikson, and J.C. Davis, “X-ray interactions: photoabsorption, scattering, transmission, and reflection at E=50-30000 eV, Z=1-92,” At. Data Nucl. Data Tables 54, 181–342 (1993). [CrossRef]

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

15. F. Arfelli, M. Assante, V. Bonvicini, A. Bravin, G. Cantatore, E. Castelli, L. D. Palma, M. D. Michiel, R. Longo, A. Olivo, S. Pani, D. Pontoni, P. Poropat, M. Prest, A. Rashevsky, G. Tromba, A. Vacchi, E. Vallazza, and F. Zanconati, “Low-dose phase contrast x-ray medical imaging,” Phys. Med. Biol. 43, 2845–2852 (1998). [CrossRef] [PubMed]

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

16. 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,” Radiology 215, 286–293 (2000).

, 17

17. M. Z. Kiss, D. E. Sayers, Z. Zhong, C. Parham, and E. D. Pisano, “Improved image contrast of calcifications in breast tissue specimens using diffraction enhanced imaging,” Phys. Med. Biol. 49, 3427–3439 (2004). [CrossRef] [PubMed]

, 18

18. E. Pisano, R. Johnston, D. Chapman, J. Geradts, M. Iacocca, C. Livasy, D. Washburn, D. Sayers, Z. Zhong, M. Kiss, and W. Thomlinson, “Human Breast Cancer Specimens: Diffraction-enhanced Imaging with Histologic Correlation-Improved Conspicuity of Lesion Detail Compared with Digital Radiography,” Radiology 214, 895–901 (2000). [PubMed]

, 19

19. T. Tanaka, C. Honda, S. Matsuo, K. Noma, H. Ohara, N. Nitta, S. Ota, K. Tsuchiya, Y. Sakashita, A. Yamada, M. Yamasaki, A. Furukawa, M. Takahashi, and K. Murata, “The first trial of phase contrast imaging for digital full-field mammography using a practical molybdenum x-ray tube,” Invest. Radiol. 40, 385–396 (2005). [CrossRef] [PubMed]

] and cartilage imaging [20

20. C. Muehleman, J. Li, Z. Zhong, J. G. Brankov, and M. N. Wernick, “Multiple-image radiography for soft tissue of the foot and ankle,” J. Anat. 208, 115–124 (2006). [CrossRef] [PubMed]

, 21

21. J. G. Brankov, M. N. Wernick, Y. Yang, J. Li, C. Muehleman, Z. Zhong, and M. A. Anastasio, “A computed tomography implementation of multiple-image radiography,” Med. Phys. 33, 278–289 (2006). [CrossRef] [PubMed]

].

Most experimental implementations of X-ray phase-contrast imaging [4

4. R. Lewis, “Medical phase contrast x-ray imaging: current status and future prospects,” Phys. Med. Biol. 49, 3573–3583 (2004). [CrossRef] [PubMed]

] that hold promise for clinical imaging are broadly characterized as being either analyzer-based or propagation-based. In analyzer-based methods [22

22. A. Bravin, “Exploiting the x-ray refraction contrast with an analyser: the state of the art,” J. Phys. D 36, A24–A29 (2003).

, 11

11. M. N. Wernick, O. Wirjadi, D. Chapman, Z. Zhong, N. P. Galatsanos, Y. Yang, J. G. Brankov, O. Oltulu, M. A. Anastasio, and C. Muehleman, “Multiple-image radiography,” Phys. Med. Biol. 48, 3875–3895 (2003). [CrossRef]

], the object is irradiated with a quasi-monochromatic, parallel-beam, X-ray wavefield, and an analyzer crystal or other X-ray diffracting element [13

13. F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Phase retrieval and differential phase-contrast imaging with low-brilliance X-ray sources,” Nature Phys. 2, 258–261 (2006). [CrossRef]

] 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

23. P. Spanne, C. Raven, I. Snigireva, and A. Snigirev, “In-line holography and phase-contrast microtomography with high energy x-rays,” Phys. Med. Biol. 44, 741–749 (1999). [CrossRef] [PubMed]

, 6

6. T. Davis, D. Gao, T. E. Gureyev, A. Stevenson, and S. Wilkins, “Phase-contrast imaging of weakly absorbing materials using hard X-rays,” Nature (London) 373, 335–338 (1996).

, 24

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

, 25

25. A. Pogany, D. Gao, and S. W. Wilkins, “Contrast and resolution in imaging with a microfocus x-ray source,” Rev. Sci. Instrum. 68, 2774–2782 (1997). [CrossRef]

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

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

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

26. M. Born and E. Wolf, Principles of Optics (Cambridge University Press, 1999).

] 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

27. D. M. Paganin, T. E. Gureyev, K. M. Pavlov, R. A. Lewis, and M. Kitchen, “Quantitative phase retrieval using coherent imaging systems with linear transfer functions,” Opt. Commun. 234, 87–105 (2004). [CrossRef]

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

28. P. Cloetens, W. Ludwig, J. Baruchel, D. Dyck, J. Landuyt, J. P. Guigay, and M. Schlenker, “Holotomography: Quantitative phase tomography with micrometer resolution using hard synchrotron radiation x rays,” Appl. Phys. lett. 75, 2912–2914 (1999). [CrossRef]

, 29

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

, 30

30. T. E. Gureyev, D. M. Paganin, G. R. Myers, Y. I. Nesterest, and S. W. Wilkins, “Phase-and-amplitude computer tomography,” Appl. Phys. Lett. 89, 034102 (2006). [CrossRef]

, 31

31. M. A. Anastasio, D. Shi, F. D. Carlo, and X. Pan, “Analytic image reconstruction in local phase-contrast tomography,” Phys. Med. Biol. 49, 121–144 (2004). [CrossRef] [PubMed]

, 21

21. J. G. Brankov, M. N. Wernick, Y. Yang, J. Li, C. Muehleman, Z. Zhong, and M. A. Anastasio, “A computed tomography implementation of multiple-image radiography,” Med. Phys. 33, 278–289 (2006). [CrossRef] [PubMed]

].

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

27. D. M. Paganin, T. E. Gureyev, K. M. Pavlov, R. A. Lewis, and M. Kitchen, “Quantitative phase retrieval using coherent imaging systems with linear transfer functions,” Opt. Commun. 234, 87–105 (2004). [CrossRef]

, 32

32. Y. I. Nesterets, T. E. Gureyev, D. Paganin, K. M. Pavlov, and S. W. Wilkins, “Quantitative diffraction-enhanced x-ray imaging of weak objects,” J. Phys. D 37, 1262–1274 (2004).

, 33

33. Y. I. Nesterets, T. E. Gureyev, K. M. Pavlov, D. M. Paganin, and S. W. Wilkins, “Combined analyser-based and propagation-based phase-contrast imaging of weak objects,” Opt. Commun. 259, 19–31 (2006). [CrossRef]

, 34

34. T. E. Gureyev, A. Pogany, D.M. Paganin, and S.W. Wilkins, “Linear algorithms for phase retrieval in the Fresnel region,” Opt. Commun. 231, 53–70 (2004). [CrossRef]

] 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

28. P. Cloetens, W. Ludwig, J. Baruchel, D. Dyck, J. Landuyt, J. P. Guigay, and M. Schlenker, “Holotomography: Quantitative phase tomography with micrometer resolution using hard synchrotron radiation x rays,” Appl. Phys. lett. 75, 2912–2914 (1999). [CrossRef]

, 35

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

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

36. Y. Huang and M. A. Anastasio, “Statistically principled use of in-line measurements in intensity diffraction tomography,” J. Opt. Soc. Am. A 24, 626–642 (2007). [CrossRef]

] 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

37. D. M. Paganin, Coherent X-Ray Optics (Oxford University Press, 2006). [CrossRef]

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

38. T. E. Gureyev, Y. I. Nesterets, D. M. Paganin, A. Pogany, and S. W. Wilkins, “Linear algorithms for phase retrieval in the Fresnel region. 2. Partially coherent illumination,” Opt. Commun. 259, 569–580 (2006). [CrossRef]

, 39

39. Y. I. Nesterets, T. E. Gureyev, and S.W. Wilkins, “Polychromaticity in the combined propagation-based/analyser-based phase-contrast imaging,” J. Phys. D 38, 4259–4271 (2005).

]. The object is characterized by its complex-valued refractive index distribution

n(r)1δ(r)+jβ(r),
(1)

μ(r)=2kβ(r),
(2)

where k2πλ is the wavenumber. Note that classic radiographic methods are sensitive only to variations in β(r), while phase-contrast methods are sensitive to variations in both δ(r) and β(r). The wavefield Uo(x,y) on the object plane, which has been transmitted through the object, is given by

Uo(x,y)=T(x,y)Ui
(3)

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

T(x,y)=M(x,y)exp[jϕ(x,y)].
(4)

A(x,y)=kdzβ(r)
(5a)
ϕ(x,y)=kdzδ(r)
(5b)

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.

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 Um(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, Um(x,y) and Uo(x,y) can be regarded as the output and input, respectively, of a linear and shift-invariant coherent imaging system [40

40. K. M. Pavlov, T. E. Gureyev, D. Paganin, Y. Nesterets, M. J. Morgan, and R. A. Lewis, “Linear systems with slowly varying transfer functions and their application to X-ray phase-contrast imaging,” J. Phys. D 37, 2746–2750 (2004).

, 41

41. J. W. Goodman, Introduction to Fourier Optics, 3rd ed. (Roberts & Company Publishers, 2004).

]. In this case,

Um(x,y)=Gm(x,y)*Uo(x,y),
(6)

where ∗ denotes a two-dimensional (2D) convolution and Gm(x,y) describes the impulse response of the system in state m. In propagation-based imaging, for example, Gm(x,y) describes the Fresnel propagator [26

26. M. Born and E. Wolf, Principles of Optics (Cambridge University Press, 1999).

], with distinct values of m corresponding to different object-to-detector distances. Alternatively, in analyzer-based imaging, Gm(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

39. Y. I. Nesterets, T. E. Gureyev, and S.W. Wilkins, “Polychromaticity in the combined propagation-based/analyser-based phase-contrast imaging,” J. Phys. D 38, 4259–4271 (2005).

], Gm(x,y) describes the net effect of both.

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

Km(x,y)1Im(x,y)Ii,
(7)

where Ii=|Ui|2 is the intensity of the incident X-ray beam. Let K̃m(u,v) denote the 2D Fourier transform (FT) of Km(x,y) defined as

K˜m(u,v)=dxdyKm(x,y)exp[j2π(ux+vy)].
(8)

In biomedical imaging applications, the conditions of |A(x,y)|≪1 and slowly varying ϕ(x,y) can often be met [10

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

, 42

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

]. Under such conditions, it can be shown that [43

43. T. E. Gureyev, G. R. Myers, Y. I. Nesterets, D. Paganin, K. M. Pavlov, and S. W. Wilkins, “Stability and locality of amplitude and phase contrast tomographies,” Proc. SPIE 6318, 63180V (2006).

]

K˜m(u,v)=2G˜ma(u,v)A˜(u,v)+2G˜mp(u,v)ϕ˜(u,v)
(9)

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

G˜ma(u,v)=12[G˜m(u,v)+G˜m*(u,v)],
(10)

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

G˜mp(u,v)=12j[G˜m(u,v)G˜m*(u,v)].
(11)

Here, G̃m(u,v), Ã(u,v), and ϕ̃(u,v) are the 2D FTs of Gm(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

43. T. E. Gureyev, G. R. Myers, Y. I. Nesterets, D. Paganin, K. M. Pavlov, and S. W. Wilkins, “Stability and locality of amplitude and phase contrast tomographies,” Proc. SPIE 6318, 63180V (2006).

] for a detailed derivation of the imaging model in Eq. (9).

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

A˜(u,v)=G˜np(u,v)K˜m(u,v)G˜mp(u,v)K˜n(u,v)2[G˜ma(u,v)G˜np(u,v)G˜na(u,v)G˜mp(u,v)]
(12a)
ϕ˜(u,v)=G˜na(u,v)K˜m(u,v)+G˜ma(u,v)K˜n(u,v)2[G˜ma(u,v)G˜np(u,v)G˜na(u,v)G˜mp(u,v)].
(12b)

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

G˜ma(u,v)G˜np(u,v)G˜na(u,v)G˜mp(u,v)=0.
(13)

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 Im(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 Im(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, (N2) 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), mn=1,2,3, will be employed to describe the estimates for the case N=3, where the subscripts denote that measurements Im(x,y) and In(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

36. Y. Huang and M. A. Anastasio, “Statistically principled use of in-line measurements in intensity diffraction tomography,” J. Opt. Soc. Am. A 24, 626–642 (2007). [CrossRef]

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

ϕ˜m,n(u,v)=αl,mϕ(u,v)ϕ˜l,m(u,v)+αl,nϕ(u,v)ϕ˜l,n(u,v)
(14)
A˜m,n(u,v)=αl,ma(u,v)A˜l,m(u,v)+αl,na(u,v)A˜l,n(u,v)
(15)

with lmn=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

αl,mϕ(u,v)=G˜na(u,v)[G˜la(u,v)G˜mp(u,v)G˜ma(u,v)G˜lp(u,v)]G˜la(u,v)[G˜ma(u,v)G˜np(u,v)G˜na(u,v)G˜mp(u,v)]
(16a)
αl,nϕ(u,v)=G˜ma(u,v)[G˜la(u,v)G˜np(u,v)G˜na(u,v)G˜lp(u,v)]G˜la(u,v)[G˜ma(u,v)G˜np(u,v)G˜na(u,v)G˜mp(u,v)],
(16b)
αl,ma(u,v)=G˜np(u,v)[G˜la(u,v)G˜mp(u,v)G˜ma(u,v)G˜lp(u,v)]G˜lp(u,v)[G˜ma(u,v)G˜np(u,v)G˜na(u,v)G˜mp(u,v)],
(16c)
αl,na(u,v)=G˜mp(u,v)[G˜la(u,v)G˜np(u,v)G˜na(u,v)G˜lp(u,v)]G˜lp(u,v)[G˜ma(u,v)G˜np(u,v)G˜na(u,v)G˜mp(u,v)],
(16d)

It can be verified that

αl,mϕ(u,v)+αl,nϕ(u,v)1
(17a)
αl,ma(u,v)+αl,na(u,v)1.
(17b)

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

ϕ˜(u,v)=ω1,2ϕ(u,v)ϕ˜1,2(u,v)+ω1,3ϕ(u,v)ϕ˜1,3(u,v)
(18)
A˜(u,v)=ω1,2a(u,v)A˜1,2(u,v)+ω1,3a(u,v)A˜1,3(u,v),
(19)

where the combination coefficients satisfy

ω1,2ϕ+ω1,3ϕ=1
(20a)
ω1,2a+ω1,3a=1.
(20b)

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:

σm,n2(u,v)Var{ϕ˜m,n(u,v)}
(21)
ρk,l;m,n(r)(u,v)+jρk,l;m,n(i)(u,v)Cov{ϕ˜k,l(u,v),ϕ˜m,n(u,v)}
(22)
Rm,n(u,v)+jIm,n(u,v)ωm,nϕ(u,v),
(23)

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

Var{ϕ˜(u,v)}=ω1,2ϕ(u,v)2σ1,22+ω1,3ϕ(u,v)2σ1,32
+2Re[ω1,2ϕ(u,v)[ω1,3ϕ(u,v)]*Cov{ϕ˜1,2(u,v),ϕ˜1,3(u,v)}]
(24)

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

Var{ϕ˜}R1,2R1,2(op)=0
(25a)
Var{ϕ˜}I1,2I1,2(op)=0,
(25b)

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

R1,2(op)(u,v)=σ1,32ρ1,2;1,3(r)σ1,22+σ1,322ρ1,2;1,3(r)
(26a)
I1,2(op)(u,v)=ρ1,2;1,3(i)σ1,22+σ1,322ρ1,2;1,3(r).
(26b)

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 R2dxdyVar{ϕ(x,y)}=R2dudνVar{ϕ˜(u,ν)} 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 Im(x,y) are acquired on distinct detector planes z=zm, where m=1,2,3. The impulse response in Eq. (6) corresponds to the Fresnel propagator

Gm(x,y)=exp[jkzm]jλzmexp[jπx2+y2λzm],
(27)

and its 2D FT yields the transfer function

G˜m(u,v)=exp[jkzmjπλzm(u2+v2)].
(28)

The corresponding reconstruction formulas [24

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

] are found by use of Eq. (12):

ϕ˜m,n(u,v)=cos(πλznf2)K˜m(u,v)+cos(πλzmf2)K˜n(u,v)Dm,n
(29)
A˜m,n(u,v)=sin(πλzmf2)K˜n(u,v)+sin(πλznf2)K˜m(u,v)Dm,n,
(30)

where f 2u 2+v 2 and

Dm,n(u,v)2sin(πλf2(zmzn)).
(31)

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

u2+v2=lλ(zmzn),
(32)

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 (uM,vM) 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

uM2+vM21λ(zmzn).
(33)

Equation (33) indicates that for a fixed detector resolution, additional poles in ϕ̃m,n(u,v) will emerge when the detector spacing (zm-zn) 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

Var{ϕ˜m,n(u,v)}=cos2(πλznf2)Var{I˜m(u,v)}+cos2(πλzmf2)Var{I˜n(u,v)}Dm,n2
(34)
Var{A˜m,n(u,v)}=sin2(πλznf2)Var{I˜m(u,v)}+sin2(πλzmf2)Var{I˜n(u,v)}Dm,n2
(35)

The joint covariances of the Fourier data are computed as

Cov{ϕ˜1,2(u,v),ϕ˜1,3(u,v)}=cos(πλz2f2)cos(πλz3f2)Var{I˜1(u,v)}D1,2D1,3
(36)
Cov{A˜1,2(u,v),A˜1,3(u,v)}=sin(πλz2f2)sin(πλz3f2)Var{I˜1(u,v)}D1,2D1,3
(37)

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

Var{ϕ˜m,n(u,v)}1Dm,n2(u,v)
(38)
Var{A˜m,n(u,v)}1Dm,n2(u,v).
(39)

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

ω1,2heur(u,v)=D1,22+α1,2ϕD2,32D1,22+D1,32+D2,32
(40a)
ω1,3heur(u,v)=D1,32+α1,3ϕD2,32D1,22+D1,32+D2,32,
(40b)

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=zm is denoted as

Im[r,s]=Im(x,y)x=rΔx,y=sΔy,
(41)

where r and s are integer-valued indices that reference detector elements, and Δxy 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

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

, 44

44. S. Lowenthal and H. Arsenault, “Image formation for coherent diffuse objects: Statistical properties,” J. Opt. Soc. Am. 60, 1478–1483 (1970). [CrossRef]

]

Im[r,s]=Im0[r,s]+nm[r,s],
(42)

where I 0 m[r,s] denotes the noiseless intensity data and the signal-dependent noise nm(r,s) has zero-mean and

Var{nm[r,s]}=(Im0[r,s])2σ2(zm),
(43)

where the real-valued quantity σ 2(zm) can depend on the detector location zm [35

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

]. We also assume

Cov{nm[r,s],nm[r,s]}=Var{nm[r,s]}δrrδssδmm,
(44)

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

45. W. D. Stanley, G. R. Dougherty, and R. Dougherty, Digital Signal Processing (Reston Publishing Company, Inc., 1984).

].

The 2D DFT [45

45. W. D. Stanley, G. R. Dougherty, and R. Dougherty, Digital Signal Processing (Reston Publishing Company, Inc., 1984).

] of Eq. (42) can be written as follows

I˜m[p,q]=Im0˜[p,q]+n˜m[p,q]
(45)

where

I˜m[p,q]=r=0N1s=0N1Im[r,s]exp[j2πN(pr+qs)]
(46)
Im0˜[p,q]=r=0N1s=0N1Im0[r,s]exp[j2πN(pr+qs)]
(47)
n˜m[p,q]=r=0N1s=0N1nm[r,s]exp[j2πN(pr+qs)],
(48)

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

Var{n˜m[p,q]}=r,r=0N1s,s=0N1exp[j2πN(p(rr)+q(ss))]Cov{nm[r,s],nm[r,s]}
=r=0N1s=0N1E{(nm[r,s])2}=r=0N1s=0N1(Im0[r,s])2σ2(zm),
(49)

I˜m(u,v)u=pΔu,v=qΔvL2N2I˜m[p,q],
(50)

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

Var{I˜m(u,v)}u=pΔu,v=qΔvL4N4Var{I˜m[p,q]}=L4N4σ2(zm)r=0N1s=0N1(Im0[r,s])2.
(51)

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

Var{ϕ˜m,n(u,v)}u=pΔu,v=qΔv
L4N4cos2(πλznf2)r=0N1s=0N1(I0[r,s;zm])2σ2(zm)+cos2(πλzmf2)r=0N1s=0N1(I0[r,s;zn])2σ2(zn)Dm,n2
(52)
Var{A˜m,n(u,v)}u=pΔu,v=qΔv
L4N4sin2(πλznf2)r=0N1s=0N1(I0[r,s;zm])2σ2(zm)+sin2(πλzmf2)r=0N1s=0N1(I0[r,s;zn])2σ2(zn)Dm,n2
(53)

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

Cov{ϕ˜1,2,ϕ˜1,3}L4N4cos(πλz2f2)cos(πλz3f2)r=0N1s=0N1(I0[r,s;z1])2σ2(z1)D1,2D1,3u=pΔu,v=qΔv
(54)
Cov{A˜1,2,A˜1,3}L4N4sin(πλz2f2)sin(πλz3f2)r=0N1s=0N1(I0[r,s;z1])2σ2(z1)D1,2D1,3u=pΔu,v=qΔv.
(55)

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:

ω1,2ϕ(u,v)=K1ϕK2ϕK1ϕ+K3ϕ2K2ϕ
(56)
ω1,2ϕ(u,v)=K1aK2aK1a+K3a2K2a
(57)

where

K1ϕ=D1,22[cos2(πλz3f2)r=0N1s=0N1(I0[r,s;z1])2σ12+cos2(πλz1f2)r=0N1s=0N1(I0[r,s;z3])2σ32]
(58a)
K2ϕ=D1,2D1,3[cos(πλz2f2)cos(πλz3f2)r=0N1s=0N1(I0[r,s;z1])2σ12]
(58b)
K3ϕ=D1,32[cos2(πλz2f2)r=0N1s=0N1(I0[r,s;z1])2σ12+cos2(πλz1f2)r=0N1s=0N1(I0[r,s;z2])2σ22],
(58c)

and

K1a=D1,22[sin2(πλz3f2)r=0N1s=0N1(I0[r,s;z1])2σ12+sin2(πλz1f2)r=0N1s=0N1(I0[r,s;z3])2σ32]
(59a)
K2a=D1,2D1,3[sin(πλz2f2)sin(πλz3f2)r=0N1s=0N1(I0[r,s;z1])2σ12]
(59b)
K3a=D1,32[sin2(πλz2f2)r=0N1s=0N1(I0[r,s;z1])2σ12+sin2(πλz1f2)r=0N1s=0N1(I0[r,s;z2])2σ22],
(59c)

where σ 2 mσ 2(zm). The corresponding forms of ωϕ 1,3(u,v) and ωa 1,3(u,v) are determined by Eq. (20).

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

ω1,2ϕ(u,v)D1,22[cos2(πλz3f2)σ12+cos2(πλz1f2)σ32]D1,2D1,3[cos(πλz2f2)cos(πλz3f2)σ12]m,n=2mn3D1,m2[cos2(πλznf2)σ12+cos2(πλz1f2)σn2]2D1D1,3[cos(πλz2f2)cos(πλz3f2)σ12].
(60)

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=zm, m=1,2,3, were considered to be behind the object. The detector contained 1024×1024 elements of dimension of 1 µm2, 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 1=19 mm, z 2=96 mm, z 3=182 mm, while the corresponding positions in Geometry ‘B’ were z 1=12 mm, z 2=38 mm, and z 3=72 mm.

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

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

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 Uo(x,y) on the object plane was computed according to Eqs. (3) and (4). Subsequently, sampled values of the wavefieldUm(x,y) on each detector plane z=zm, m=1,2,3, was computed by use of Eq. (6) with Gm(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 Im[r, s] were then computed as the square of the wavefield modulus on each detector plane.

Noisy measurement data Im[r,s] were computed according to Eq. (42). The noise process nm[r,s] was described by a uncorrelated Gaussian distribution whose variance was determined by Eq. (43) with σ 2(zm)=0.05%. Ensembles of 1000 noisy realizations of Im[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 σ 2(zm) 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.
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).
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.
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).

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.

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

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

table-icon
View This 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.

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

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.

Appendix: Generalization to ≥3 measurement states

Consider a measurement geometry that consists of M measurement states 1, 2, 3, ⋯, M. Estimates of the desired object properties can be computed by use of any measurement-state pair; thus there will be ((M2)) distinct estimates available for the M-state system. Consequently, the final estimate of ϕ̃(u,v), for example, can be obtained from a weighted summation of all the available estimates. Let ϕ̃m,n(u,v) denote the phase estimate computed by use of the intensity data acquired on measurement pair (m,n). A final unbiased estimate that possesses a potentially reduced variance can be written as

ϕ˜(u,v)=m=1M1n=m+1Mω̂m,nϕ(u,v)ϕ̂m,n(u,v),
(A1)

where ω̂ϕ m,n(u,v) are generally the complex-valued coefficients that satisfy

m=1M1n=m+1Mω̂m,nϕ(u,v)=1.
(A2)

As described in Section 3, however, any estimate ϕ̃m,n(u,v) is a linearly combination of any other two estimates ϕ̃1,m(u,v) and ϕ̃1,n(u,v), where m, n are integer-valued indices. Therefore, Eq. (A1) can be simplified to

ϕ˜(u,v)=m=2Mω1,mϕ(u,v)ϕ˜1,m(u,v).
(A3)

The variance of the final estimate ϕ̃(u,v) can be obtained readily as

Var{ϕ˜(u,v)}=l=2Mω1,lϕ(u,v)2Var{ϕ˜1,l(u,v)}
+2Re[m=2M1n=m+1Mω1,mϕ(u,v)ω1,nϕ*(u,v)Cov{ϕ˜1,m(u,v)ϕ˜1,n(u,v)}],
(A4)

where l=2,3, ⋯,M in the first summation term, and m=2,3, ⋯,M-1, n=3,4, ⋯,M with n>m in the second summation term.

We introduce the following notation:

σ1,m2(u,v)Var{ϕ˜1,m[u,v]},
(A5)
ρ1,m;1,n(r)(u,v)+jρ1,m;1,n(i)(u,v)Cov{ϕ˜1,m[u,v],ϕ˜1,n[u,v]},
(A6)
R1,m(u,v)+jI1,m(u,v)ω1,mϕ(u,v).
(A7)

Thus the variance of the final Fourier estimate can be obtained,

Var{ϕ˜}=l=2M(R1,l2+I1,l2)σ1,l2
+2{m=2M1n=m+1M[ρ1,m;1,n(r)(R1,mR1,n+I1,mI1,n)ρ1,m;1,n(i)(R1,nI1,mR1,mI1,n)]}.
(A8)

In order to reach an optimal final estimate, variance needs to be specified such that,

Var{ϕ˜}R1,mR1,m(op)=0
Var{ϕ˜}I1,mI1,m(op)=0,
(A9)

where R (op) 1,m and I (op) 1,m specify the value of the optimal combination coefficient for ωϕ 1,m with m=2,3, ⋯,M. Analytic formula for determination of the combination coefficients can be derived as follows. From knowledge of Eq. (A8) and taking use of Eq. (A9), a (2M-4)×(2M-4) system is formed, which can be expressed by a matrix equation.

Hx=b
(A10)

where

H=(h11h12h1,2M4h2M4,1h2M4,2h2M4,2M4),
(A11)
x=(R1,2(op)I1,2(op)R1,3(op)I1,3(op)R1,M1(op)I1,M1(op))T,
(A12)

and

b=(b1b2b2M4)T,
(A13)

where the m-th equation in the system specifies the partial derivative of Var{ϕ̃} with respect to the m-th component of Eq. (A12). Assuming H is nonsingular, inversion of Eq. (A10) is accomplished readily:

xm(u,v)=det(Hm)det(H)m=1,2,3,,2M4
(A14)

where H m is the matrix obtained by replacing the m-th column of H with b. If we specify the coefficient ωϕ 1,m(u,v) to be real-valued, the term I 1,m(u,v) will vanish, yielding a (M-2)×(M-2) system.

Analogous to Section 4.3, the heuristically determined combination coefficient ω heur 1,m for M- detector system can be computed readily as:

ω1,mheur(u,v)=D1,m2+α1,mn=2nmMDm,n2l=1M1m=l+1MDl,m2
(A15)

Acknowledgment

This research was supported in part by National Institutes of Health research grant EB004507 and a National Science Foundation CAREER Award 0546113.

References and links

1.

A. Krol, A. Ikhlef, J.-C. Kieffer, D. Bassano, C. C. Chamberlain, Z. Jiang, H. Pepin, and S. C. Prasad, “Laser-based microfocused X-ray source for mammography: feasibility study,” Med. Phys. 24, 725–732 (1997). [CrossRef] [PubMed]

2.

R. Waynant, “Toward practical coherent x-ray sources: Potential medical applications,” IEEE J. Quantum Electron. 6, 1465–1469 (2000). [CrossRef]

3.

H. Yamada, “Novel x-ray source based on a tabletop synchrotron and its unique features,” Nucl. Instrum. Methods Phys. Res. B 199, 509–516 (2003).

4.

R. Lewis, “Medical phase contrast x-ray imaging: current status and future prospects,” Phys. Med. Biol. 49, 3573–3583 (2004). [CrossRef] [PubMed]

5.

W. Thomlinson, P. Suortti, and D. Chapman, “Recent advances in synchrotron radiation medical research,” Nucl. Instrum. Methods Phys. Res. A 543, 288–296 (2005).

6.

T. Davis, D. Gao, T. E. Gureyev, A. Stevenson, and S. Wilkins, “Phase-contrast imaging of weakly absorbing materials using hard X-rays,” Nature (London) 373, 335–338 (1996).

7.

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

8.

D. Chapman, W. Thomlinson, R. E. Johnston, D. Washburn, E. Pisano, N. Gmur, Z. Zhong, R. Menk, F. Arfelli, and D. Sayers, “Diffraction enhanced x-ray imaging,” Phys. Med. Biol. 42, 2015–2025 (1997). [CrossRef] [PubMed]

9.

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]

10.

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]

11.

M. N. Wernick, O. Wirjadi, D. Chapman, Z. Zhong, N. P. Galatsanos, Y. Yang, J. G. Brankov, O. Oltulu, M. A. Anastasio, and C. Muehleman, “Multiple-image radiography,” Phys. Med. Biol. 48, 3875–3895 (2003). [CrossRef]

12.

E. F. Donnelly, R. R. Price, and D. R. Pickens, “Characterization of the phase-contrast radiography edge-enhancement effect in a cabinet x-ray system,” Med. Phys. 30, 2292–2296 (2003). [CrossRef] [PubMed]

13.

F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Phase retrieval and differential phase-contrast imaging with low-brilliance X-ray sources,” Nature Phys. 2, 258–261 (2006). [CrossRef]

14.

B.L. Henke, E.M. Gullikson, and J.C. Davis, “X-ray interactions: photoabsorption, scattering, transmission, and reflection at E=50-30000 eV, Z=1-92,” At. Data Nucl. Data Tables 54, 181–342 (1993). [CrossRef]

15.

F. Arfelli, M. Assante, V. Bonvicini, A. Bravin, G. Cantatore, E. Castelli, L. D. Palma, M. D. Michiel, R. Longo, A. Olivo, S. Pani, D. Pontoni, P. Poropat, M. Prest, A. Rashevsky, G. Tromba, A. Vacchi, E. Vallazza, and F. Zanconati, “Low-dose phase contrast x-ray medical imaging,” Phys. Med. Biol. 43, 2845–2852 (1998). [CrossRef] [PubMed]

16.

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,” Radiology 215, 286–293 (2000).

17.

M. Z. Kiss, D. E. Sayers, Z. Zhong, C. Parham, and E. D. Pisano, “Improved image contrast of calcifications in breast tissue specimens using diffraction enhanced imaging,” Phys. Med. Biol. 49, 3427–3439 (2004). [CrossRef] [PubMed]

18.

E. Pisano, R. Johnston, D. Chapman, J. Geradts, M. Iacocca, C. Livasy, D. Washburn, D. Sayers, Z. Zhong, M. Kiss, and W. Thomlinson, “Human Breast Cancer Specimens: Diffraction-enhanced Imaging with Histologic Correlation-Improved Conspicuity of Lesion Detail Compared with Digital Radiography,” Radiology 214, 895–901 (2000). [PubMed]

19.

T. Tanaka, C. Honda, S. Matsuo, K. Noma, H. Ohara, N. Nitta, S. Ota, K. Tsuchiya, Y. Sakashita, A. Yamada, M. Yamasaki, A. Furukawa, M. Takahashi, and K. Murata, “The first trial of phase contrast imaging for digital full-field mammography using a practical molybdenum x-ray tube,” Invest. Radiol. 40, 385–396 (2005). [CrossRef] [PubMed]

20.

C. Muehleman, J. Li, Z. Zhong, J. G. Brankov, and M. N. Wernick, “Multiple-image radiography for soft tissue of the foot and ankle,” J. Anat. 208, 115–124 (2006). [CrossRef] [PubMed]

21.

J. G. Brankov, M. N. Wernick, Y. Yang, J. Li, C. Muehleman, Z. Zhong, and M. A. Anastasio, “A computed tomography implementation of multiple-image radiography,” Med. Phys. 33, 278–289 (2006). [CrossRef] [PubMed]

22.

A. Bravin, “Exploiting the x-ray refraction contrast with an analyser: the state of the art,” J. Phys. D 36, A24–A29 (2003).

23.

P. Spanne, C. Raven, I. Snigireva, and A. Snigirev, “In-line holography and phase-contrast microtomography with high energy x-rays,” Phys. Med. Biol. 44, 741–749 (1999). [CrossRef] [PubMed]

24.

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

25.

A. Pogany, D. Gao, and S. W. Wilkins, “Contrast and resolution in imaging with a microfocus x-ray source,” Rev. Sci. Instrum. 68, 2774–2782 (1997). [CrossRef]

26.

M. Born and E. Wolf, Principles of Optics (Cambridge University Press, 1999).

27.

D. M. Paganin, T. E. Gureyev, K. M. Pavlov, R. A. Lewis, and M. Kitchen, “Quantitative phase retrieval using coherent imaging systems with linear transfer functions,” Opt. Commun. 234, 87–105 (2004). [CrossRef]

28.

P. Cloetens, W. Ludwig, J. Baruchel, D. Dyck, J. Landuyt, J. P. Guigay, and M. Schlenker, “Holotomography: Quantitative phase tomography with micrometer resolution using hard synchrotron radiation x rays,” Appl. Phys. lett. 75, 2912–2914 (1999). [CrossRef]

29.

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

30.

T. E. Gureyev, D. M. Paganin, G. R. Myers, Y. I. Nesterest, and S. W. Wilkins, “Phase-and-amplitude computer tomography,” Appl. Phys. Lett. 89, 034102 (2006). [CrossRef]

31.

M. A. Anastasio, D. Shi, F. D. Carlo, and X. Pan, “Analytic image reconstruction in local phase-contrast tomography,” Phys. Med. Biol. 49, 121–144 (2004). [CrossRef] [PubMed]

32.

Y. I. Nesterets, T. E. Gureyev, D. Paganin, K. M. Pavlov, and S. W. Wilkins, “Quantitative diffraction-enhanced x-ray imaging of weak objects,” J. Phys. D 37, 1262–1274 (2004).

33.

Y. I. Nesterets, T. E. Gureyev, K. M. Pavlov, D. M. Paganin, and S. W. Wilkins, “Combined analyser-based and propagation-based phase-contrast imaging of weak objects,” Opt. Commun. 259, 19–31 (2006). [CrossRef]

34.

T. E. Gureyev, A. Pogany, D.M. Paganin, and S.W. Wilkins, “Linear algorithms for phase retrieval in the Fresnel region,” Opt. Commun. 231, 53–70 (2004). [CrossRef]

35.

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]

36.

Y. Huang and M. A. Anastasio, “Statistically principled use of in-line measurements in intensity diffraction tomography,” J. Opt. Soc. Am. A 24, 626–642 (2007). [CrossRef]

37.

D. M. Paganin, Coherent X-Ray Optics (Oxford University Press, 2006). [CrossRef]

38.

T. E. Gureyev, Y. I. Nesterets, D. M. Paganin, A. Pogany, and S. W. Wilkins, “Linear algorithms for phase retrieval in the Fresnel region. 2. Partially coherent illumination,” Opt. Commun. 259, 569–580 (2006). [CrossRef]

39.

Y. I. Nesterets, T. E. Gureyev, and S.W. Wilkins, “Polychromaticity in the combined propagation-based/analyser-based phase-contrast imaging,” J. Phys. D 38, 4259–4271 (2005).

40.

K. M. Pavlov, T. E. Gureyev, D. Paganin, Y. Nesterets, M. J. Morgan, and R. A. Lewis, “Linear systems with slowly varying transfer functions and their application to X-ray phase-contrast imaging,” J. Phys. D 37, 2746–2750 (2004).

41.

J. W. Goodman, Introduction to Fourier Optics, 3rd ed. (Roberts & Company Publishers, 2004).

42.

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

43.

T. E. Gureyev, G. R. Myers, Y. I. Nesterets, D. Paganin, K. M. Pavlov, and S. W. Wilkins, “Stability and locality of amplitude and phase contrast tomographies,” Proc. SPIE 6318, 63180V (2006).

44.

S. Lowenthal and H. Arsenault, “Image formation for coherent diffuse objects: Statistical properties,” J. Opt. Soc. Am. 60, 1478–1483 (1970). [CrossRef]

45.

W. D. Stanley, G. R. Dougherty, and R. Dougherty, Digital Signal Processing (Reston Publishing Company, Inc., 1984).

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

ToC Category:
Imaging Systems

History
Original Manuscript: February 5, 2007
Revised Manuscript: May 7, 2007
Manuscript Accepted: May 20, 2007
Published: July 25, 2007

Virtual Issues
Vol. 2, Iss. 9 Virtual Journal for Biomedical Optics

Citation
Cheng-Ying Chou, Yin Huang, Daxin Shi, and Mark A. Anastasio, "Image reconstruction in quantitative X-ray phase-contrast imaging employing multiple measurements," Opt. Express 15, 10002-10025 (2007)
http://www.opticsinfobase.org/oe/abstract.cfm?URI=oe-15-16-10002


Sort:  Year  |  Journal  |  Reset  

References

  1. A. Krol, A. Ikhlef, J.-C. Kieffer, D. Bassano, C. C. Chamberlain, Z. Jiang, H. Pepin, and S. C. Prasad, "Laserbased microfocused X-ray source for mammography: feasibility study," Med. Phys. 24, 725-732 (1997). [CrossRef] [PubMed]
  2. R. Waynant, "Toward practical coherent x-ray sources: Potential medical applications," IEEE J. Quantum Electron. 6, 1465-1469 (2000). [CrossRef]
  3. H. Yamada, "Novel x-ray source based on a tabletop synchrotron and its unique features," Nucl. Instrum.Methods Phys. Res. B 199, 509-516 (2003).
  4. R. Lewis, "Medical phase contrast x-ray imaging: current status and future prospects," Phys. Med. Biol. 49, 3573-3583 (2004). [CrossRef] [PubMed]
  5. W. Thomlinson, P. Suortti, and D. Chapman, "Recent advances in synchrotron radiation medical research," Nucl. Instrum. Methods Phys. Res. A 543, 288-296 (2005).
  6. T. Davis, D. Gao, T. E. Gureyev, A. Stevenson, and S. Wilkins, "Phase-contrast imaging of weakly absorbing materials using hard X-rays," Nature (London) 373, 335-338 (1996).
  7. K. A. Nugent, T. E. Gureyev, D. Cookson, D. Paganin, and Z. Barnea, "Quantitative phase imaging using hard x-rays," Phys. Rev. Lett. 77, 2961-2964 (1996). [CrossRef] [PubMed]
  8. D. Chapman, W. Thomlinson, R. E. Johnston, D. Washburn, E. Pisano, N. Gmur, Z. Zhong, R. Menk, F. Arfelli, and D. Sayers, "Diffraction enhanced x-ray imaging," Phys. Med. Biol. 42, 2015-2025 (1997). [CrossRef] [PubMed]
  9. 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]
  10. 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]
  11. M. N. Wernick, O. Wirjadi, D. Chapman, Z. Zhong, N. P. Galatsanos, Y. Yang, J. G. Brankov, O. Oltulu, M. A. Anastasio, and C. Muehleman, "Multiple-image radiography," Phys. Med. Biol. 48, 3875-3895 (2003). [CrossRef]
  12. E. F. Donnelly, R. R. Price, and D. R. Pickens, "Characterization of the phase-contrast radiography edgeenhancement effect in a cabinet x-ray system," Med. Phys. 30, 2292-2296 (2003). [CrossRef] [PubMed]
  13. F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, "Phase retrieval and differential phase-contrast imaging with low-brilliance X-ray sources," Nature Phys. 2, 258-261 (2006). [CrossRef]
  14. B.L. Henke, E.M. Gullikson, and J.C. Davis, "X-ray interactions: photoabsorption, scattering, transmission, and reflection at E=50-30000 eV, Z=1-92," At. Data Nucl. Data Tables 54, 181-342 (1993). [CrossRef]
  15. F. Arfelli, M. Assante, V. Bonvicini, A. Bravin, G. Cantatore, E. Castelli, L. D. Palma, M. D. Michiel, R. Longo, A. Olivo, S. Pani, D. Pontoni, P. Poropat, M. Prest, A. Rashevsky, G. Tromba, A. Vacchi, E. Vallazza, and F. Zanconati, "Low-dose phase contrast x-ray medical imaging," Phys. Med. Biol. 43, 2845-2852 (1998). [CrossRef] [PubMed]
  16. 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," Radiology 215, 286-293 (2000).
  17. M. Z. Kiss, D. E. Sayers, Z. Zhong, C. Parham, and E. D. Pisano, "Improved image contrast of calcifications in breast tissue specimens using diffraction enhanced imaging," Phys. Med. Biol. 49, 3427-3439 (2004). [CrossRef] [PubMed]
  18. E. Pisano, R. Johnston, D. Chapman, J. Geradts, M. Iacocca, C. Livasy, D. Washburn, D. Sayers, Z. Zhong, M. Kiss, and W. Thomlinson, "Human Breast Cancer Specimens: Diffraction-enhanced Imaging with Histologic Correlation-Improved Conspicuity of Lesion Detail Compared with Digital Radiography," Radiology 214, 895- 901 (2000). [PubMed]
  19. T. Tanaka, C. Honda, S. Matsuo, K. Noma, H. Ohara, N. Nitta, S. Ota, K. Tsuchiya, Y. Sakashita, A. Yamada, M. Yamasaki, A. Furukawa, M. Takahashi, and K. Murata, "The first trial of phase contrast imaging for digital full-field mammography using a practical molybdenum x-ray tube," Invest. Radiol. 40, 385-396 (2005). [CrossRef] [PubMed]
  20. C. Muehleman, J. Li, Z. Zhong, J. G. Brankov, and M. N. Wernick, "Multiple-image radiography for soft tissue of the foot and ankle," J. Anat. 208, 115-124 (2006). [CrossRef] [PubMed]
  21. J. G. Brankov, M. N. Wernick, Y. Yang, J. Li, C. Muehleman, Z. Zhong, and M. A. Anastasio, "A computed tomography implementation of multiple-image radiography," Med. Phys. 33, 278-289 (2006). [CrossRef] [PubMed]
  22. A. Bravin, "Exploiting the x-ray refraction contrast with an analyser: the state of the art," J. Phys. D 36, A24-A29 (2003).
  23. P. Spanne, C. Raven, I. Snigireva, and A. Snigirev, "In-line holography and phase-contrast microtomography with high energy x-rays," Phys. Med. Biol. 44, 741-749 (1999). [CrossRef] [PubMed]
  24. 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).
  25. A. Pogany, D. Gao, and S. W. Wilkins, "Contrast and resolution in imaging with a microfocus x-ray source," Rev. Sci. Instrum. 68, 2774-2782 (1997). [CrossRef]
  26. M. Born and E. Wolf, Principles of Optics (Cambridge University Press, 1999).
  27. D. M. Paganin, T. E. Gureyev, K. M. Pavlov, R. A. Lewis, and M. Kitchen, "Quantitative phase retrieval using coherent imaging systems with linear transfer functions," Opt. Commun. 234, 87-105 (2004). [CrossRef]
  28. P. Cloetens, W. Ludwig, J. Baruchel, D. Dyck, J. Landuyt, J. P. Guigay, and M. Schlenker, "Holotomography: Quantitative phase tomography with micrometer resolution using hard synchrotron radiation x rays," Appl. Phys. lett. 75, 2912-2914 (1999). [CrossRef]
  29. A. V. Bronnikov, "Theory of quantitative phase-contrast computed tomography," J. Opt. Soc. Am. A 19, 472-480 (2002). [CrossRef]
  30. T. E. Gureyev, D. M. Paganin, G. R. Myers, Y. I. Nesterest, and S. W. Wilkins, "Phase-and-amplitude computer tomography," Appl. Phys. Lett. 89, 034102 (2006). [CrossRef]
  31. M. A. Anastasio, D. Shi, F. D. Carlo, and X. Pan, "Analytic image reconstruction in local phase-contrast tomography," Phys. Med. Biol. 49, 121-144 (2004). [CrossRef] [PubMed]
  32. Y. I. Nesterets, T. E. Gureyev, D. Paganin, K. M. Pavlov, and S. W. Wilkins, "Quantitative diffraction-enhanced x-ray imaging of weak objects," J. Phys. D 37, 1262-1274 (2004).
  33. Y. I. Nesterets, T. E. Gureyev, K. M. Pavlov, D. M. Paganin, and S. W. Wilkins, "Combined analyser-based and propagation-based phase-contrast imaging of weak objects," Opt. Commun. 259, 19-31 (2006). [CrossRef]
  34. T. E. Gureyev, A. Pogany, D.M. Paganin, and S.W. Wilkins, "Linear algorithms for phase retrieval in the Fresnel region," Opt. Commun. 231, 53-70 (2004). [CrossRef]
  35. 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]
  36. Y. Huang and M. A. Anastasio, "Statistically principled use of in-line measurements in intensity diffraction tomography," J. Opt. Soc. Am. A 24, 626-642 (2007). [CrossRef]
  37. D. M. Paganin, Coherent X-Ray Optics (Oxford University Press, 2006). [CrossRef]
  38. T. E. Gureyev, Y. I. Nesterets, D. M. Paganin, A. Pogany, and S. W. Wilkins, "Linear algorithms for phase retrieval in the Fresnel region. 2. Partially coherent illumination," Opt. Commun. 259, 569-580 (2006). [CrossRef]
  39. Y. I. Nesterets, T. E. Gureyev, and S.W. Wilkins, "Polychromaticity in the combined propagation-based/analyserbased phase-contrast imaging," J. Phys. D 38, 4259-4271 (2005).
  40. K. M. Pavlov, T. E. Gureyev, D. Paganin, Y. Nesterets, M. J. Morgan, and R. A. Lewis, "Linear systems with slowly varying transfer functions and their application to X-ray phase-contrast imaging," J. Phys. D 37, 2746- 2750 (2004).
  41. J. W. Goodman, Introduction to Fourier Optics, 3rd ed. (Roberts & Company Publishers, 2004).
  42. J.-P. Guigay, "Fourier transform analysis of Fresnel diffraction patterns and in-line holograms," Optik 49, 121- 125 (1977).
  43. T. E. Gureyev, G. R. Myers, Y. I. Nesterets, D. Paganin, K. M. Pavlov, and S. W. Wilkins, "Stability and locality of amplitude and phase contrast tomographies," Proc. SPIE 6318, 63180V (2006).
  44. S. Lowenthal and H. Arsenault, "Image formation for coherent diffuse objects: Statistical properties," J. Opt. Soc. Am. 60, 1478-1483 (1970). [CrossRef]
  45. W. D. Stanley, G. R. Dougherty, and R. Dougherty, Digital Signal Processing (Reston Publishing Company, Inc., 1984).

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