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

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 p hase-contrast imaging is the presence of isolated Fourier domain singular ities, which can greatly amplify the noise levels in the estimated Fourie r domain and lead to noisy and/or distorted images in spatial domain. In t his article, we develop a statistically optimal reconstruction method tha t employs multiple (> 2) measurement states to mitigate the noise amplification ef fects due to singularities in the reconstruction formula. Computer-si mulation studies are carried out to quantitatively and systematically investig ate the developed method, within the context of propagation-based X-ray phas e-contrast imaging. The reconstructed images are shown to possess dram atic lly reduced noise levels and greatly enhanced imaging contrast . © 2007 Optical Society of America OCIS codes:(110.7440) X-ray imaging; (100.5070) Phase retrieval; (170 .3010) Image reconstruction techniques. References and links 1. A. Krol, A. Ikhlef, J.-C. Kieffer, D. Bassano, C. C. Chamber lain, Z. Jiang, H. Pepin, and S. C. Prasad, “Laserbased microfocused X-ray source for mammography: feasibility s tudy,” Med. Phys.24, 725–732 (1997). 2. R. Waynant, “Toward practical coherent x-ray sources: Po tential medical applications,” IEEE J. Quantum Electron.6, 1465–1469 (2000). 3. H. Yamada, “Novel x-ray source based on a tabletop synchrot ron and its unique features,” Nucl. Instrum. Methods Phys. Res. B199, 509–516 (2003). 4. R. Lewis, “Medical phase contrast x-ray imaging: current s ta us and future prospects,” Phys. Med. Biol. 49, 3573–3583 (2004). 5. W. Thomlinson, P. Suortti, and D. Chapman, “Recent advances i 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. Wilkin s, “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. B arnea, “Quantitative phase imaging using hard x-rays,” Phys. Rev. Lett. 77, 2961–2964 (1996). 8. D. Chapman, W. Thomlinson, R. E. Johnston, D. Washburn, E. Pi sano, N. Gmur, Z. Zhong, R. Menk, F. Arfelli, and D. Sayers, “Diffraction enhanced x-ray imaging,” Phys. M ed. Biol.42, 2015–2025 (1997). 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). 10. X. Wu and H. Liu, “Clinical implementation of X-ray phase-c ontrast imaging: Theoretical foundations and design considerations,” Med. Phys. 30, 2169–2179 (2003). 11. M. N. Wernick, O. Wirjadi, D. Chapman, Z. Zhong, N. P. Galat sanos, Y. Yang, J. G. Brankov, O. Oltulu, M. A. Anastasio, and C. Muehleman, “Multiple-image radiography,” Phys. Med. Biol.48, 3875–3895 (2003). #79616 $15.00 USD Received 5 Feb 2007; revised 7 May 2007; accepted 20 May 2007; published 25 Jul 2007 (C) 2007 OSA 6 August 2007 / Vol. 15, No. 16 / OPTICS EXPRESS 10002 12. E. F. Donnelly, R. R. Price, and D. R. Pickens, “Character ization of the phase-contrast radiography edgeenhancement effect in a cabinet x-ray system,” Med. Phys. 30, 2292–2296 (2003). 13. F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Phase ret rieval and differential phase-contrast imaging with low-brilliance X-ray sources,” Nature Phys. 2, 258–261 (2006). 14. B.L. Henke, E.M. Gullikson, and J.C. Davis, “X-ray inter actions: photoabsorption, scattering, transmission, and reflection at E=50-30000 eV, Z=1-92,” At. Data Nucl. Data Tab les54, 181–342 (1993). 15. F. Arfelli, M. Assante, V. Bonvicini, A. Bravin, G. Canta tore, E. Castelli, L. D. Palma, M. D. Michiel, R. Longo, A. Olivo, S. Pani, D. Pontoni, P. Poropat, M. Prest, A. Rashev ky, G. Tromba, A. Vacchi, E. Vallazza, and F. Zanconati, “Low-dose phase contrast x-ray medical imaging ,” Phys. Med. Biol.43, 2845–2852 (1998). 16. F. Arfelli, V. Bonvicini, A. Bravin, G. Cantatore, E. Cas telli, 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 Sy nchrotron Radiation: Phase-Detection Techniques,” Radiology215, 286–293 (2000). 17. M. Z. Kiss, D. E. Sayers, Z. Zhong, C. Parham, and E. D. Pisan o, “Improved image contrast of calcifications in breast tissue specimens using diffraction enhanced imaging, ” Phys. Med. Biol.49, 3427–3439 (2004). 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: Dif fraction-enhanced Imaging with Histologic Correlation-Improved Conspicuity of Lesion Detail Compared with Digital Radiography,” Radiology214, 895– 901 (2000). 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 fir st t ial of phase contrast imaging for digital full-field mammography using a practical molybdenum x-ray tube, ” Invest. Radiol.40, 385–396 (2005). 20. C. Muehleman, J. Li, Z. Zhong, J. G. Brankov, and M. N. Werni ck, “Multiple-image radiography for soft tissue of the foot and ankle,” J. Anat. 208, 115–124 (2006). 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). 22. A. Bravin, “Exploiting the x-ray refraction contrast wi th 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-l ine holography and phase-contrast microtomography with high energy x-rays,” Phys. Med. Biol. 44, 741–749 (1999). 24. P. Cloetens, “Contribution to Phase Contrast Imaging, Re construction and Tomography with Hard Synchrotron Radiation: Principles, Implementation and Applications,” P h.D. thesis, Vrije Universiteit Brussel (1999). 25. A. Pogany, D. Gao, and S. W. Wilkins, “Contrast and resolut i n in imaging with a microfocus x-ray source,” Rev. Sci. Instrum. 68, 2774–2782 (1997). 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,” Op t. Commun.234, 87–105 (2004). 28. P. Cloetens, W. Ludwig, J. Baruchel, D. Dyck, J. Landuyt, J. P. Guigay, and M. Schlenker, “Holotomography: Quantitative phase tomography with micrometer resolution usi g hard synchrotron radiation x rays,” Appl. Phys. lett. 75, 2912–2914 (1999). 29. A. V. Bronnikov, “Theory of quantitative phase-contras computed tomography,” J. Opt. Soc. Am. A 19, 472–480 (2002). 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). 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). 32. Y. I. Nesterets, T. E. Gureyev, D. Paganin, K. M. Pavlov, a nd 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 , a d S. W. Wilkins, “Combined analyser-based and propagation-based phase-contrast imaging of weak objects, ” Opt. Commun.259, 19–31 (2006). 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). 35. D. Paganin, A. Barty, P. J. Mcmahon, and K. A. Nugent, “Quan tit tive phase-amplitude microscopy III. The effects of noise,” J. Microsc. 214, 51–61 (2003). 36. Y. Huang and M. A. Anastasio, “Statistically principled use of in-line measurements in intensity diffraction tomography,” J. Opt. Soc. Am. A24, 626–642 (2007). 37. D. M. Paganin, Coherent X-Ray Optics (Oxford University Press, 2006). 38. T. E. Gureyev, Y. I. Nesterets, D. M. Paganin, A. Pogany, a d S. W. Wilkins, “Linear algorithms for phase retrieval in the Fresnel region. 2. Partially coherent illu mination,” Opt. Commun. 259, 569–580 (2006). 39. Y. I. Nesterets, T. E. Gureyev, and S. W. Wilkins, “Polych romaticity 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 #79616 $15.00 USD Received 5 Feb 2007; revised 7 May 2007; accepted 20 May 2007; published 25 Jul 2007 (C) 2007 OSA 6 August 2007 / Vol. 15, No. 16 / OPTICS EXPRESS 10003 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 di ffraction 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 cohere nt diffuse objects: Statistical properties,” J. Opt. Soc. Am.60, 1478–1483 (1970). 45. W. D. Stanley, G. R. Dougherty, and R. Dougherty, Digital Signal Processing(Reston Publishing Company, Inc., 1984).


Introduction
Due to the advent and accessibility of synchrotron and laboratory X-ray sources [1,2,3] that possess high degrees of coherence, new imaging methods, broadly characterized as X-ray phase-contrast imaging methods [4,5,6,7,8,9,10,11,12,13] are being developed that have dramatic advantages over conventional radiographic and tomographic X-ray imaging systems.Unlike conventional radiographic methods, phase-contrast imaging methods exploit a contrast mechanism based on differences in the complex-valued X-ray refractive index values of tissue.At diagnostic X-ray energies, variations in the real component of the refractive index of tissues are several orders of magnitude larger than are variations in the imaginary component, or equivalently, the X-ray attenuation coefficient [14].Consequently, phase-contrast imaging may permit the visualization and quantification of tissues that have identical, or very similar, X-ray absorption properties.Additionally, unlike absorption-contrast, the phase-contrast mechanism persists at relatively high energies, which permits low-dose imaging [15].Experimental studies of X-ray phase-contrast imaging have revealed significantly enhanced contrast and tissue discriminability over conventional radiographic methods in applications including cancer detection [16,17,18,19] and cartilage imaging [20,21].
Most experimental implementations of X-ray phase-contrast imaging [4] that hold promise for clinical imaging are broadly characterized as being either analyzer-based or propagationbased.In analyzer-based methods [22,11], the object is irradiated with a quasi-monochromatic, parallel-beam, X-ray wavefield, and an analyzer crystal or other X-ray diffracting element [13] is present between the object and detector.The analyzer crystal diffracts the components of the transmitted wavefield traveling at or near the crystal's Bragg angle, thereby rejecting all wavefield components traveling outside a narrow angular range.This angular filtering of the transmitted wavefield produces a phase-contrast enhancement in the recorded intensity image.Propagation-based methods [23,6,24,25], also known as in-line methods, do not require the use of diffractive X-ray elements, and employ a classic in-line holographic measurement geometry that is essentially similar to magnification radiography [9].The object is irradiated with a plane-wave or paraxial X-ray wavefield possessing a sufficient degree of spatial coherence, and the intensity of the transmitted wavefield is recorded by a detector placed at some non-zero distance behind the object.In this case, the measured intensity distribution represents an in-line (Gabor) hologram [26] that contains coded information about the refractive properties of the object.The phase shifts that are introduced into the probing wavefield by the refractive index variations within the object are transferred into intensity variations in the final measurement by the process of free-space (Fresnel) wave propagation between the object and detector.
Quantitative X-ray phase-contrast imaging methods seek to reconstruct separate images that depict the object's projected absorption and real-valued refractive index distributions, which reflect two distinct and complementary intrinsic object properties.To achieve this, the measured intensity data must generally be recorded at two or more distinct "states" of the imaging system [27].For example, in analyzer-based methods, different measurement states correspond to distinct orientations of the analyzer crystal, while in propagation-based methods they could correspond to distinct object-to-detector distances.Quantitative phase-contrast imaging methods are computed-imaging methods, and require use of reconstruction algorithms for image formation.When implemented in tomographic mode, an estimate of the three-dimensional (3D) complex refractive index distribution can be reconstructed by use of phase-contrast tomography methods [28,29,30,31,21].
Under certain assumptions regarding the transmitted wavefield, phase-contrast imaging systems can be modeled as linear, shift-invariant, imaging systems.Assuming measurements are acquired at two distinct states of the imaging system, Fourier-based reconstruction formulas can be derived [27,32,33,34] from knowledge of the imaging system's optical transfer function.A common feature of these reconstruction formulas, however, is the presence of isolated Fourier domain singularities.In practice, the Fourier components of the projected object properties residing near the singularities can contain greatly amplified noise levels, resulting in noisy and distorted images.Because the locations of the singularities are determined by the chosen two measurement states, their effects can be mitigated by acquiring intensity measurements at multiple (> 2) measurement states [28,35].For example, different pairs of intensity measurements can be utilized to reconstruct different Fourier components of the object properties, which provides the opportunity to potentially avoid the singularities.However, such simple strategies do not exploit fully the statistically complementary information contained in the available measurement data, and can result in images with suboptimal statistical properties.
In this work, we propose and investigate a methodology for image reconstruction in quantitative phase-contrast imaging when multiple (> 2) intensity measurements are available.Linear estimators are proposed that combine the available intensity measurements to produce statistical estimates of the projected object properties whose Fourier components possess optimally reduced variances.This general strategy is inspired by a recent study of intensity diffraction tomography [36] by our group.Explicit forms of the reconstruction formulas are derived for propagation-based phase-contrast imaging, where a general noise model and finite-sampling effects are considered.Computer-simulation studies are conducted to demonstrate the efficacy of the proposed method.

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

Interaction of X-ray wavefield with object
As depicted in Fig. 1, consider that an object is irradiated by a monochromatic X-ray wavefield U i with wavelength λ , which is traveling along the positive z-axis.The effects of imperfect wavefield coherence will not be considered, but can be addressed as in [38,39].The object is characterized by its complex-valued refractive index distribution where r = (x, y, z), and δ ( r) and β ( r) quantify the X-ray refractive and absorption properties of the object.The quantity β ( r) is related to the linear X-ray attenuation coefficient µ( r) as where k ≡ 2π λ 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 U o (x, y) on the object plane, which has been transmitted through the object, is given by where T (x, y) is the transmission function that can be expressed generally as The amplitude modulus M(x, y) = exp[−A(x, y)] and the phase shift φ (x, y) describe how the amplitude and wavefront (i.e., phase), respectively, of the probing wavefield are perturbed by the presence of the object.They are related to imaginary and real components, respectively, of n( r) as where the line-integrals are computed over the support of the object.
Fig. 1.A schematic of a generic X-ray phase-contrast imaging system.

Linear shift-invariant X-ray phase-contrast imaging systems
Consider again the generic X-ray phase-contrast imaging system shown in Fig. 1.Let U m (x, y) denote the transmitted wavefield on a detector plane of constant z, which is downstream from the object plane.The integer-valued subscript m is employed to denote the state of the imaging system.For many analyzer-and propagation-based phase-contrast imaging systems, U m (x, y) and U o (x, y) can be regarded as the output and input, respectively, of a linear and shift-invariant coherent imaging system [40,41].In this case, where * denotes a two-dimensional (2D) convolution and G m (x, y) describes the impulse response of the system in state m.In propagation-based imaging, for example, G m (x, y) describes the Fresnel propagator [26], with distinct values of m corresponding to different object-todetector distances.Alternatively, in analyzer-based imaging, G m (x, y) describes the coherent impulse response of the analyzer crystal or other diffractive element(s) employed, with distinct values of m corresponding to different orientations of the analyzer.In hybrid systems [39], G m (x, y) describes the net effect of both.
On the detector plane, the intensity I m (x, y) = |U m (x, y)| 2 is recorded, which represents a radiograph with mixed absorption-and phase-contrast.From knowledge of the measured intensity I m (x, y), a modified data function can be defined as where is the intensity of the incident X-ray beam.Let K m (u, v) denote the 2D Fourier transform (FT) of K m (x, y) defined as In biomedical imaging applications, the conditions of |A(x, y)| << 1 and slowly varying φ (x, y) can often be met [10,42].Under such conditions, it can be shown that [43] where G a m (u, v) is the amplitude transfer function (ATF): and G p m (u, v) is the phase transfer function (PTF): Here, G m (u, v), A(u, v), and φ (u, v) are the 2D FTs of G m (x, y), A(x, y), and φ (x, y), respectively, and G * m (•, •) denotes the complex conjugate of G m (•, •).The interested reader is referred to Ref. [43] for a detailed derivation of the imaging model in Eq. (9).
Equation (9) relates the measured intensity I m (x, y), or equivalently K m (x, y), to the 2D FTs of the sought-after quantities A(x, y) and φ (x, y).If an additional measurement I n (x, y) is obtained when the imaging system is in state n = m, A(u, v) and φ (u, v) can be determined algebraically as Equations (12a) and (b) represent reconstruction formulas for quantitative phase-contrast imaging.From the determined A(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 When Eq. ( 12) is applied to noisy, or otherwise inconsistent, measurement data, the Fourier components (u, v) of A(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.

Variance reduction in quantitative X-ray phase-contrast imaging
We consider that intensity measurements I m (x, y) are acquired at three distinct states m = 1, 2, 3 of the system.The results below are generalized to the case of an arbitrary number of measurements in the Appendix.The intensity data function I m (x, y) is interpreted as a stochastic process, which reflects that the measurements are contaminated by stochastic errors such as detector noise.Our goal is to exploit the statistically complementary information in the available measurements to reduce the variance of the estimated A(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, N 2 estimates of A(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 A m,n (u, v) and φ m,n (u, v), m = n = 1, 2, 3, will be employed to describe the estimates for the case N = 3, where the subscripts denote that measurements I m (x, y) and I n (x, y) were employed.Because the locations of poles in Eq. ( 12) depend on the choice of measurement states, the components of A 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 A m,n (u, v) or φ m,n (u, v) in a way that attempts to cancel the poles in each two-state estimate [36], to form final estimates A(u, v) or φ (u, v), respectively, that possess reduced variances for all (u, v).However, it should be noted that only two of the three available estimates are independent in the sense that with l = m = n = 1, 2, 3.In other words, any estimate can be expressed as a linear combination of the remaining two.The coefficients in Eqs. ( 14) and ( 15) are frequency-dependent and given by Accordingly, final estimates of A(u, v) and φ (u, v) that exploit statistically complementary information in the three intensity measurements can be formed as where the combination coefficients satisfy Equations ( 18) and ( 19) represent linear estimators for producing unbiased estimates of φ (u, v) and A(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 A m,n (u, v), respectively.Moreover, as described next, they can be designed to optimally reduce the variance of φ (u, v) and A(u, v).We consider first the problem of producing estimates φ (u, v) having reduced variances.The following notation will be employed: where 'Var' and 'Cov' denote the variance and covariance, respectively of a random process.
where the superscript * denotes the complex conjugate.In order to minimize Var{ φ (u, v)}, we need that where R 1,2 and I 1,2 are the real and imaginary components of ω φ 1,2 , respectively.The solution of these equations yields The optimal choice of ω φ 1,3 (u, v) is subsequently determined by Eq. (20a).The specification of the combination coefficients ω a 1,2 (u, v) and ω a 1,3 (u, v) that optimally reduce the variance of A(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 R 2 dxdyVar{φ (x, y)} = R 2 dudvVar{ φ (u, v)} 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).

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.

Two-state reconstruction formulas
In propagation-based imaging, the different system states can correspond to different objectto-detector distances.We consider that three intensity measurements I m (x, y) are acquired on distinct detector planes z = z m , where m = 1, 2, 3.The impulse response in Eq. ( 6) corresponds to the Fresnel propagator and its 2D FT yields the transfer function The corresponding reconstruction formulas [24] are found by use of Eq. ( 12): where As before, the indices m, n satisfy m = 1, 2, n = 2, 3 with n > m.
Note that D m,n = 0 is equivalent to Eq. ( 13), and specifies the locations of poles in the reconstruction formulas.Equation ( 29) contains poles corresponding to frequencies (u, v) that satisfy where l is an integer.One such pole is located at zero-frequency u = v = 0, indicating that the low-frequency components of φ (u, v) will contain highly amplified noise levels.The existence of additional poles, away from the origin of the Fourier space, depends on the detector resolution and the detector pair spacing.Let (u M , v M ) denote the maximum spatial frequencies recorded by the detector.Additional poles in the reconstruction formulas described by Eqs. ( 29) and (30) will be present when Equation (33) indicates that for a fixed detector resolution, additional poles in φ m,n (u, v) will emerge when the detector spacing (z m − z n ) is sufficiently large.Likewise, this discussion of poles also applies to A m,n (u, v), with the exception that the pole at u = v = 0 is not present due to a cancellation.

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 A(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 A(u, v).
From Eqs. ( 29) and ( 30), one finds readily that The joint covariances of the Fourier data are computed as It should be noted that Eqs. ( 36) and ( 37) are real-valued, and therefore the imaginary components of the optimal combination coefficients will vanish [see Eq. (26b)].

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 A m,n (u, v) can still be effectively mitigated by suitable heuristic specification of the combination coefficients.From Eqs. ( 34)- (35), near the locations of poles, we find that Heuristic combination coefficients ω heur m,n (u, v) should have a (u, v)-dependence that is inversely proportional to that indicated in Eqs. ( 38)- (39).Moreover, at locations of poles we should have ω heur m,n (u, v) ≡ 0. Based on these requirements, the ω heur m,n (u, v) for use in estimating φ (u, v) via Eq.( 18) can be chosen as where α φ 1,2 and α φ 1,3 are defined in Eq. ( 16).When estimating A(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).

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.

Noise model and finite sampling considerations
The discretely sampled intensity data on a measurement plane z = z m is denoted as where r and s are integer-valued indices that reference detector elements, and ∆x = ∆y denotes the element dimension in a square detector array of dimension L × L. Equation ( 41) assumes idealized (Dirac delta) sampling, namely, the averaging effects of sampling aperture are not considered.However, the analysis follows can be generalized to address such effects.The square bracket, '[•]', is introduced to represent the functions whose arguments are discretely sampled.
We assume the noise model satisfies [35,44] where I 0 m [r, s] denotes the noiseless intensity data and the signal-dependent noise n m (r, s) has zero-mean and where the real-valued quantity σ 2 (z m ) can depend on the detector location z m [35].We also assume where δ mn denotes the Kronecker delta function.

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 I m (u, v) is required.To compute this from knowledge of discretely sampled data, the continuous 2D FT will be approximated by use of the discrete Fourier transform (DFT) [45].The 2D DFT [45] of Eq. ( 42) can be written as follows where with p, q denoting the integer-valued Fourier indices that are conjugate to [r, s], and N specifying the number of detector elements in each dimension of the square 2D detector.The variance of n m [p, q] is computed as where 'E{•}' denotes the statistical expectation operator.The continuous and discrete FTs of the intensity data are related as [45] where ∆u = ∆v = 1 L are the frequency domain sampling intervals along u and v axes.By use of Eqs. ( 49) and (50), we find that Finally, using Eq. ( 51) with the results of Section 4.2, we arrive at Var{ φ m,n (u, v)} u=p∆u,v=q∆v where, as before, m = 1, 2, n = 2, 3 with n > m, and
If ∑ N−1 r=0 ∑ N−1 s=0 I 0 [r, s; z m ] 2 does not vary significantly as a function of m (i.e., the detector location z m ), Eq. ( 56) can be expressed in the somewhat simplified form: (60)

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

Numerical phantom and in-line measurement geometry
The in-line measurement geometry shown in Fig. 2 was assumed.A monochromatic X-ray plane-wave with wavelength λ = 1 × 10 −10 m, propagated along the z-axis and irradiated an object.Three detector planes located at z = z m , m = 1, 2, 3, were considered to be behind the object.The detector contained 1024×1024 elements of dimension of 1 µm 2 , 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.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.

Measurement data and simulation studies
The simulated intensity measurements were computed as follows.From knowledge of φ (x, y) and A(x, y), the transmitted wavefield U o (x, y) on the object plane was computed according to Eqs. ( 3) and (4).Subsequently, sampled values of the wavefield U m (x, y) on each detector plane z = z m , m = 1, 2, 3, was computed by use of Eq. ( 6) with G m (x, y) specified by the Fresnel propagator in Eq. ( 27).The convolution in Eq. ( 6) was computed by use of the 2D fast Fourier transform (FFT) algorithm.The intensity data I m [r, s] were then computed as the square of the wavefield modulus on each detector plane.
Noisy measurement data I m [r, s] were computed according to Eq. ( 42).The noise process n m [r, s] was described by a uncorrelated Gaussian distribution whose variance was determined by Eq. ( 43) with σ 2 (z m ) = 0.05%.Ensembles of 1000 noisy realizations of I m [r, s] were com- puted for m = 1, 2, 3.

Image reconstruction
Estimates φ m,n (u, v) and A 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 A 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 A 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 A(u, v) that possess optimally reduced variances.The required combination coefficients were computed according to Eqs. ( 56) and (57).Because σ 2 (z m ) was fixed at a constant value, it can be verified that, in this special case, the optimal combination coefficients given in Eqs. ( 56) and ( 57) are identical to the heuristic ones defined by Eq. ( 40).Corresponding estimates of φ (x, y) and A(x, y) were computed by application of the 2D inverse FFT algorithm.The variances of the reconstructed object properties in both the Fourier and spatial domains were estimated empirically.

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{ A(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{ A(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{ A(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{ A(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.
Variance data for the images contained in Fig. 8 are shown in Fig. 9. Figures 9(a) and (b) display empirical variance estimates of the reconstructed Fourier components.Each subfigure contains profiles corresponding to the Fourier variance estimated by use of detector pair (1,2) (dashed curve), detector pair (1,3) (dashdotted curve), detector (2,3) (dotted curve), and the optimal ones reconstructed by use of Eqs. ( 18) or ( 19) (solid curve).As expected, the optimal estimates possess variances that are lower than any of the two-state estimates, for all frequency components.Subfigures (c) and (d) display empirical estimates of the corresponding images in the spatial domain.The labeling of the profiles is the same as described above.The variances of the optimally estimated images are seen to be lower than the two-state reconstructions, for all values of (x, y).This reflects that estimates having reduced Fourier variances will generally have reduced variances in the spatial domain.
The corresponding results obtained with Geometry 'B' are shown in Figs. 10 and 11.From Fig. 10, we find again that the optimal estimates of A(x, y) and φ (x, y) contain obviously reduced noise levels as compared to the two-state reconstructions.The quantitative variance data in Fig. 11 confirms that the optimal estimates of A(x, y) and φ (x, y) contain lower variances than the two-state reconstructions.
In practice, the detector locations cannot be determined with perfect precision.The uncertainty in the detector positions will introduce inconsistencies between the measured data and assumed imaging model.In order to investigate the effects of this, images were reconstructed from the noisy intensity measurements corresponding to Geometry 'A', where imperfect knowledge of the detector positions was assumed.The assumed detector positions ẑm were related to the true locations z m by ẑm = z m + ε(z m ), with m = 1, 2, 3, where ε(z m ) represents the positioning errors.The three sets of errors ε(z m ) contained in Table 1 were considered.The corresponding estimates of A(x, y) and φ (x, y) reconstructed by the optimal estimation method are shown in Fig. 12. Subfigures (a)-(c) contain the optimal estimates of A(x, y) for detector position error levels (1), (2), and (3), respectively.The corresponding estimates of φ (x, y) are contained in subfigures (d)-(f).Despite the geometry uncertainties, the reconstructed A(x, y) and φ (x, y) closely resemble the corresponding images in Fig. 8(d) and (h), which were reconstructed from the same noisy data but assuming perfect knowledge of the imaging geometry.The mean square error (MSE) values of the estimates of A(x, y) displayed in Fig. 12 crete (Cartesian grid) representation of the true functions A(x, y) and φ (x, y), and do not reflect an average over the distribution of the measurement noise.These observations suggest that our multi-plane reconstruction method may be robust to geometry errors under certain conditions.However, a detailed investigation of the effects of geometry errors on the proposed image reconstruction methods remains an important topic for future research.

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 a methodology 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.
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 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 ( M 2 ) 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

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

Table 1 .
Error levels in the detector positions in Geometry 'A'.