## Abstract

Quantitative in-line X-ray phase-contrast imaging methods seek to reconstruct separate images that depict an object’s projected absorption and refractive properties. An understanding of the statistical properties of the reconstructed images can facilitate the identification of optimal imaging parameters for specific diagnostic tasks. However, the statistical properties of quantitative X-ray phase-contrast imaging remain largely unexplored. In this work, we derive analytic expressions that describe the second-order statistics of the reconstructed absorption and phase images. Concepts from statistical decision theory are applied to demonstrate how the statistical properties of images corresponding to distinct imaging geometries can influence signal detectability.

© 2009 Optical Society of America

## 1. Introduction

In-line X-ray phase-contrast imaging methods [1–7] are being actively developed and possess significant advantages over conventional X-ray imaging methods. Because they exploit a contrast mechanism based on differences in an object’s X-ray refractive index distribution, they can permit the visualization of features that possess identical, or very similar, X-ray absorption properties. Moreover, they can operate effectively at relatively high X-ray energies, thereby facilitating low-dose imaging [8]. An attractive feature of in-line, or propagation-based, phase-contrast imaging is that it possesses a relatively simple implementation that is essentially similar to magnification radiography, and attains phase-contrast enhancement via free space wave propagation between the object and detector.

By use of specialized data acquisition protocols and phase retrieval formulas [4, 9–12], it is possible to perform quantitative phase-contrast imaging in which separate images depicting the projected absorption and refractive properties of an object are reconstructed. Phase retrieval formulas based on the transport of intensity equation (TIE) [1,4,13,14] and the contrast transfer function (CTF) [15, 16] have been developed and investigated. Due to the underlying physics, the phase retrieval formulas unavoidably possess a pole at the origin in Fourier space, resulting in an amplification of low-frequency noise [17] in the reconstructed phase image.

Physical measures of image quality such as noise variance, signal contrast, and spatial resolution metrics, have been applied to optimize X-ray phase-contrast imaging systems [17–22]. In the broader context of modern image science, however, there is growing consensus that image quality should be specified and measured by the average performance of an observer on a specific diagnostic task [23–25]. Task-based measures of image quality, also referred to as objective measures, are fundamentally distinct from physical measures in that they are inherently statistical and are grounded in statistical decision theory. Computation of task-based measures of image quality requires knowledge of the statistical properties of an image, which corresponds to a realization of a stochastic process. It is known [24, 26, 27], for example, that the correlation structure of the noise, and not simply the noise variance and signal contrast, can strongly influence the ability of a human or numerical observer to detect a specified signal in a noisy image. There remains an important need to conduct a quantitative analysis of the statistical properties of images in X-ray phase-contrast imaging, which would facilitate the use of task-based image quality measures to further develop and optimize this emerging method for specific applications.

In this work, we derive analytic expressions that describe the second-order statistics of the reconstructed absorption and phase images in their continuous and discrete forms. Due to the zero-frequency pole in the phase retrieval formula, the reconstructed phase image is demonstrated to possess strong spatial correlations that are not present in the reconstructed absorption image or a conventional X-ray radiograph. The strength of these image correlations is shown to depend on the imaging geometry. Knowledge of the images’ covariance matrices and concepts from statistical decision theory are subsequently applied to demonstrate how task-based image quality, as measured by signal detectability in a signal-known-exactly/background-known-exactly (SKE/BKE) detection task, can be influenced by the imaging geometry. In this initial study, we do not address the effects of partial beam coherence, but our statistical analyses can be readily generalized to do so.

## 2. Background: in-line X-ray phase-contrast imaging

Below, the salient features of in-line X-ray phase-contrast imaging are reviewed briefly. We refer the reader to the monograph by Paganin [1] and the literature cited therein for a comprehensive treatment of phase-contrast image formation.

#### 2.1. Measurement model

We assume the canonical in-line imaging geometry depicted in Fig. 1. An object that is characterized by its X-ray refractive index distribution *n*(*r*⃗)≡1-*δ*(*r*⃗)+*jβ*(*r*⃗), *j*≡√-1, is irradiated by monochromatic and time-harmonic X-ray plane wave*U _{i}* with wavelength

*λ*. The irradiating wave propagates in the direction of the positive

*z*-axis. Immediately behind the object on the object plane

*z*=0, the transmitted wavefield

*U*(

_{t}*x,y*) can be described as

where *ϕ* (*x,y*) and *A*(*x,y*) are the phase and amplitude perturbations introduced by the object. Assuming the validity of the projection approximation [1],* ϕ* (*x,y*) and *A*(*x,y*) are related to *n*(*r*⃗) as

and

where $k\equiv \genfrac{}{}{0.1ex}{}{2\pi}{\lambda}$ denotes the wavenumber and *L* denotes the paths of X-ray beams, which are approximated to be parallel with the *z*-axis.

Let *U _{m}*(

*x,y*) denote the transmitted wavefield on a detector plane

*z*=

*z*, which is down-stream from the object plane. Note that the integer-valued subscript

_{m}*m*is employed to index the object-to-detector distance

*z*. The wavefield intensity data recorded on this detector plane are given by

_{m}where *h _{m}*(

*x,y*) denotes the associated Fresnel propagator kernel and *

_{2}indicates a two-dimensional (2D) convolution operation. From knowledge of the measured intensity

*I*(

_{m}*x,y*), we define the data function

where *I _{i}*=|

*U*|

_{i}^{2}is the intensity of the incident X-ray beam.

The goal of quantitative in-line phase-contrast imaging is to reconstruct estimates of *ϕ*(*x,y*) and *A*(*x,y*) from knowledge of two or more intensity measurements *I _{m}*(

*x,y*) acquired at distinct object-to-detector distances. More generally, the intensity measurements could be acquired on a fixed detector plane and an alternative degree of freedom in the imaging system could be varied [28–30] to acquire the necessary measurement data.

#### 2.2. Image reconstruction based on the contrast transfer function (CTF)

Let

and

denote the 2D Fourier transforms (FTs) of the sought after quantities *ϕ* (*x,y*) and *A*(*x,y*). We assume that the two intensity measurements *I _{m}*(

*x,y*) and

*I*(

_{n}*x,y*) are acquired on the detector planes

*z*=

*z*and

_{m}*z*=

*z*, respectively, and their 2D FTs will be denoted as

_{n}*Ĩ*(

_{m}*u,v*) and

*Ĩ*(

_{n}*u,v*). Similarly, the 2D FT of the data function

*K*(

_{m}*x,y*) will be denoted as

*K*̃

_{m}(

*u,v*).

Based on the CTF formulation of phase-contrast image formation [6, 12], the Fourier components of *ϕ *(*x,y*) and *A*(*x,y*) can be determined as

and

where *f _{2}*≡

*u*+

_{2}*v*and the subscripts (

_{2}*m,n*) have been added to

*ϕ*̃(

*u,v*) and

*A*̃(

*u,v*) to denote that they have been determined by the intensity measurements

*I*(

_{m}*x,y*) and

*I*(

_{n}*x,y*). Subsequently, estimates of

*ϕ*(

_{m,n}*x,y*) and Am,n(

*x,y*) can be obtained via the inverse 2D FT. In many applications, including medical imaging [2], the detector pixel size and imaging geometry establish the near field conditions sin[

*πλ f*

^{2}(

*z*)]≈

_{n}-z_{m}*πλ f*

^{2}(

*z*) and cos[

_{n}-z_{m}*πλ f*

^{2}

*z*]≈cos[

_{m}*πλ f*

^{2}

*z*]≈1. For simplicity we will assume such conditions, under which Eqs. (8) and (9) can be simplified as

_{n}and

However, that the statistical analyses below can be immediately extended to the case where Eqs. (8) and (9) must be considered.

We note that Eq. (10) has a pole at zero frequency (*u*=0,*v*=0), while Eq. (11) does not. As investigated next, this will cause the reconstructed estimates of *ϕ _{m,n}*(

*x,y*) and

*A*(

_{m,n}*x,y*) to possess markedly different second-order statistical properties that will affect signal detection tasks.

## 3. Second-order statistics of image estimates: Continuous case

Below we derive closed-form expressions for the covariance of the estimated phase and absorption images. Our initial analysis is based on a continuous formulation of the problem, but finite sampling effects are subsequently addressed in Section 4. Without a change in notation, the intensity data function *I _{m}*(

*x,y*), and hence the data function

*K*(

_{m}*x,y*) and reconstructed images, will be interpreted as stochastic processes, reflecting that in practice they are inherently random quantities. We will let Cov{

*A,B*} and Var{

*A*} denote the covariance and variance of the stochastic processes

*A*and

*B*[23, 31].

#### 3.1. Autocovariance of the Fourier estimates

**Phase image:** By use of the definition of the autocovariance [31] and the phase retrieval formula in Eq. (10), it can be verified that

where *f*′^{2}≡*u*′^{2}+*v*′^{2} and

$$=\int {\int}_{\infty}\mathbf{d}x\mathbf{d}y\int {\int}_{\infty}\mathbf{d}x\prime \mathbf{d}y\prime \mathrm{exp}\left[-j2\pi \left[\left(u-u\prime \right)x+\left(v-v\prime \right)y\right]\right]\mathrm{Cov}\left\{{K}_{m}(x,y),{K}_{m}(x\prime ,y\prime )\right\}.$$

For the case of uncorrelated noise,

where *δ* (·) and *δ _{m,m}*′ denote the 1D Dirac and Kronecker delta functions, and therefore Eq. (12) can be expressed as

$$=\genfrac{}{}{0.1ex}{}{1}{4{\pi}^{2}{\lambda}^{2}{\left({z}_{n}-{z}_{m}\right)}^{2}{f}^{2}{f\prime}^{2}}$$

$$\times \int {\int}_{\infty}\mathbf{d}x\mathbf{d}y\mathrm{exp}\left[-j2\pi \left[\left(u-u\prime \right)x+(v-v\prime )y\right]\right]\left[\mathrm{Var}\left\{{K}_{m}(x,y)\right\}+\mathrm{Var}\left\{{K}_{n}(x,y)\right\}\right].$$

An important feature of Eq. (12) or (15) is the Fourier space pole located at *f*=*f*′=0. This pole indicates that the low frequency components of *ϕ*
_{m,n}(*x,y*) will be highly correlated. As demonstrated in our numerical studies, this will cause the noise in the reconstructed estimates of the phase images to possess a lumpy appearance.

**Absorption image:** In an analogous way, by use of Eq. (11) the Fourier autocovariance of the absorption image can be determined as

$$\times \left[{z}_{n}^{2}\mathrm{Cov}\left\{{\tilde{K}}_{m}(u,v),{\tilde{K}}_{m}(u\text{'},v\text{'})\right\}+{z}_{m}^{2}\mathrm{Cov}\left\{{\tilde{K}}_{n}(u,v),{\tilde{K}}_{n}(u\text{'},v\text{'})\right\}\right].$$

In the case of uncorrelated noise as described by Eq. (14), Eq. (16) can be expressed as

$$=\genfrac{}{}{0.1ex}{}{1}{{4\left({z}_{n}-{z}_{m}\right)}^{2}}\int {\int}_{\infty}\mathbf{d}x\mathbf{d}y\mathrm{exp}\left[-j2\pi \left[\left(u-u\prime \right)x+\left(v-v\prime \right)y\right]\right]\left[{z}_{n}^{2}\mathrm{Var}\left\{{K}_{m}(x,y)\right\}+{z}_{m}^{2}\mathrm{Var}\left\{{K}_{n}(x,y)\right\}\right].$$

Unlike the Fourier autocovariance for the phase image in Eq. (15), Eq. (17) does not possess a pole at the origin of Fourier space (or anywhere else). Accordingly, we expect that the reconstructed estimate of the absorption image will not contain strong low-frequency correlations and will possess a noise texture that is markedly different than that of the phase image.

#### 3.2. Autocovariance of the image estimates

From knowledge of its Fourier autocovariance functions, the autocovariance of the phase and absorption images can be readily determined as

$$\times \int {\int}_{\infty}\mathbf{d}u\prime \mathbf{d}v\prime \mathrm{exp}\left[-j2\pi (u\prime x\prime +v\prime y\prime )\right]\mathrm{Cov}\{{\tilde{\varphi}}_{m,n}(u,v),{\tilde{\varphi}}_{m,n}(u\prime ,v\prime )\},$$

and

$$\times \int {\int}_{\infty}\mathbf{d}u\prime \mathbf{d}v\prime \mathrm{exp}\left[-j2\pi (u\prime x\prime +v\prime y\prime )\right]\mathrm{Cov}\{{\tilde{A}}_{m,n}(u,v),{\tilde{A}}_{m,n}(u\prime ,v\prime )\}.$$

Below, we derive explicit forms of Eqs. (18) and (19) with consideration of finite sampling effects, and utilize the resulting expressions in a study of signal detectability.

## 4. Second-order statistics of image estimates: Discrete case

#### 4.1. Finite sampling considerations and discrete noise model

When a digital detector is employed, the measured intensity data should be mathematically described as an ordered collection of numbers (i.e., a vector) rather than a function of a continuous variable. We will denote the sampled intensity data on a measurement plane *z*=*z _{m}* as

where *r* and *s* are integer-valued detector-element indices and Δ*x*=Δ*y* denotes the element dimension in a square detector array of dimension *L*×*L*. Equation (20) 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. Here and elsewhere, a square bracket, ‘[·]’, is employed to describe discretely sampled quantities, in order to distinguish them from functions of continuous variables.

We assume the measured intensity data satisfy [17, 32]

where *I _{m}*[

*r, s*] and

*I*

^{0}

*[*

_{m}*r, s*] denote the noise-contaminated and noiseless measurement data, respectively. The noise term

*n*[

_{m}*r, s*] is assumed to have a zero mean and could be signal-dependent. We also assume the noise satisfies

where *σ*
^{2}[*r, s; z _{m}*]≡Var{

*n*[

_{m}*r, s*]}. This implies that the autocovariance of the sampled data function

*K*[

_{m}*r, s*]=

*I*[

_{m}*r, s*]/

*I*-1 satisfies

_{i}where we have regarded the intensity *I _{i}* of the incident beam as a deterministic quantity.

#### 4.2. Second-order statistics

In order to compute the second-order statistics described in Section 3, knowledge of the autocovariance of *K*̃* _{m}*(

*u,v*) is required. To compute this in the case of of discretely sampled measurement data, the continuous 2D FT will be approximated by use of the discrete Fourier transform (DFT) [33].

The 2D DFT of Eq. (21) can be written as follows

where

and

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. Analogously, the 2D DFT of the sampled data function *K _{m}*[

*r, s*] will be denoted as

*K*̃

*[*

_{m}*p,q*].

The continuous and discrete FTs of the data function are related as [33]

where $\mathbf{\Delta}u=\mathbf{\Delta}v=\genfrac{}{}{0.1ex}{}{1}{L}$ are the frequency domain sampling intervals along *u *and *v* axes. In order for this approximation to be useful, we consider the effects of aliasing to be negligible. A systematic investigation of the role of pre-sampling blur on image statistics and signal detectability in phase-contrast image formation remains an important topic for future studies [34]. By use of Eqs. (13)–(14), (23) and (28), the sampled values of Fourier autocovariance of the data function can be determined as

$$=\genfrac{}{}{0.1ex}{}{{L}^{4}}{{N}^{4}}\underset{r=0}{\overset{N-1}{\Sigma}}\underset{s=0}{\overset{N-1}{\Sigma}}\mathrm{exp}\left[-j\genfrac{}{}{0.1ex}{}{2\pi}{N}\left(r\left(p-p\prime \right)+s\left(q-q\prime \right)\right)\right]\mathrm{Var}\left\{{K}_{m}[r,s]\right\}$$

$$=\genfrac{}{}{0.1ex}{}{{L}^{4}}{{N}^{4}{I}_{i}^{2}}\underset{r=0}{\overset{N-1}{\Sigma}}\underset{s=0}{\overset{N-1}{\Sigma}}\mathrm{exp}\left[-j\genfrac{}{}{0.1ex}{}{2\pi}{N}\left(r\left(p-p\prime \right)+s\left(q-q\prime \right)\right)\right]{\sigma}^{2}[r,s,{z}_{m}].$$

Consequently, by use of Eqs. (12) and (28)–(29), the Fourier autocovariance of the sampled phase images can be determined as,

$$=\genfrac{}{}{0.1ex}{}{{L}^{4}}{{N}^{4}}\genfrac{}{}{0.1ex}{}{\mathrm{Cov}\left\{{\tilde{K}}_{m}[p,q],{\tilde{K}}_{m}[p\prime ,q\prime ]\right\}+\mathrm{Cov}\left\{{\tilde{K}}_{n}[p,q],{\tilde{K}}_{n}[p\prime ,q\prime ]\right\}}{4{\pi}^{2}{\lambda}^{2}{\left({z}_{n}-{z}_{m}\right)}^{2}{f}_{d}^{2}{f\prime}_{d}^{2}},$$

where ${f}_{d}^{2}\equiv \left({p}^{2}+{q}^{2}\right)\genfrac{}{}{0.1ex}{}{1}{{L}^{2}}$and ${f\prime}_{d}^{2}\equiv \left({p\prime}^{2}+{q}^{2}\right)\genfrac{}{}{0.1ex}{}{1}{{L}^{2}}$. By use of the inverse DFT, the autocovariance of the phase image can be computed as

$$=\genfrac{}{}{0.1ex}{}{1}{{N}^{4}}\underset{p=0}{\overset{N-1}{\mathbf{\Sigma}}}\underset{q=0}{\overset{N-1}{\mathbf{\Sigma}}}\mathrm{exp}\left[j\genfrac{}{}{0.1ex}{}{2\pi}{N}\left(\mathrm{pr}+\mathrm{qs}\right)\right]$$

$$\times \underset{p\text{'}=0}{\overset{N-1}{\mathbf{\Sigma}}}\underset{q\text{'}=0}{\overset{N-1}{\mathbf{\Sigma}}}\mathrm{exp}\left[-j\genfrac{}{}{0.1ex}{}{2\pi}{N}\left(p\prime r\prime +q\prime s\prime \right)\right]\mathrm{Cov}\{{\tilde{\varphi}}_{m,n}[p,q],{\tilde{\varphi}}_{m,n}[p\prime ,q\prime ]\}.$$

The same strategy can be applied to determine the Fourier and image autocovariances in the the absorption case, where we find that

$$=\genfrac{}{}{0.1ex}{}{{L}^{4}}{{N}^{4}}\genfrac{}{}{0.1ex}{}{1}{{4\left({z}_{n}-{z}_{m}\right)}^{2}}\left[{z}_{n}^{2}\mathrm{Cov}\left\{{\tilde{K}}_{m}[p,q],{\tilde{K}}_{m}[p\prime ,q\prime ]\right\}+{z}_{m}^{2}\mathrm{Cov}\left\{{\tilde{K}}_{n}[p,q],{\tilde{K}}_{n}[p\prime ,q\prime ]\right\}\right],$$

and

$$=\genfrac{}{}{0.1ex}{}{1}{{N}^{4}}\underset{p=0}{\overset{N-1}{\mathbf{\Sigma}}}\underset{q=0}{\overset{N-1}{\mathbf{\Sigma}}}\mathrm{exp}\left[j\genfrac{}{}{0.1ex}{}{2\pi}{N}\left(\mathrm{pr}+\mathrm{qs}\right)\right]$$

$$\times \underset{p\prime =0}{\overset{N-1}{\mathbf{\Sigma}}}\underset{q\prime =0}{\overset{N-1}{\mathbf{\Sigma}}}\mathrm{exp}\left[-j\genfrac{}{}{0.1ex}{}{2\pi}{N}\left(p\prime r\prime +q\prime s\prime \right)\right]\mathrm{Cov}\{{\tilde{A}}_{m,n}[p,q],{\tilde{A}}_{m,n}[p\prime ,q\prime ]\}.$$

Below, the expressions for the autocovariances of the phase and absorption images described by Eqs. (31) and (33) are employed with concepts from statistical decision theory to demonstrate how the image noise properties in quantitative phase-contrast imaging can affect signal detectability.

## 5. Task-based image quality: A signal detection study

In modern image science, it is well accepted [23,35] that image quality should be defined in the context of the intended use of the images. For example, in medical imaging applications, image quality is associated with a specific diagnostic task such as tumor detection. Image quality is affected by measurement noise, object variability, the physical characteristics of the imaging system, as well as the skill of the observer to use the images to make a correct diagnosis. Task-based measures of image quality measure the average performance of the observer computed over an ensemble of noisy images. In practice, this observer is often a human, but it is often too time-consuming and/or expensive to utilize humans in the repeated assessments of image quality that are required when optimizing an imaging system. To circumvent this difficulty, numerical observer models [25] are commonly employed as surrogate human observers. In our studies, we employed a well-studied numerical observer called the Hotelling observer [23, 36] to detect signals in noisy simulated absorption and phase images corresponding to different imaging geometries. For a review of numerical observers and applications of statistical decision theory in image science, we refer the reader to the monograph by Barrett and Myers [23].

At this point, it is useful to introduce some additional notation. Let the random vector $\mathbf{g}\in {\mathbb{R}}^{{N}^{2}}$ denote a lexicographically ordered representation of the image *ϕ _{m,n}*[

*r, s*] or

*A*[

_{m,n}*r, s*], where

*N*is the dimension of the digital image. The matrix

**K**will denote the image autocovariance E{(

_{g}**g**-

**g**̄)(

**g**-

**g**̄)†}, where † is the conjugate transpose (adjoint) operation, E{·} denotes the statistical expectation operator associated with the ensemble of noisy images, and

**g**̄≡E{

**g**} is the mean image. Note that the elements of

**K**are specified by Eqs. (31) or (33).

_{g}#### 5.1. Signal detection task and the Hotelling discriminant

We consider a simple signal-known-exactly/background-known-exactly (SKE/BKE) detection task. In this task, a known signal (described below) and known background object are considered. The task of the observer is to make a binary decision as to whether the signal is present or absent in each collection of noisy images presented to it. We employed the numerical Hotelling observer to make this decision. The Hotelling observer provides optimal signal detection performance for the class of observers that are restricted to performing only linear operations on the data [23]. Let the images with the signal present and signal absent be denoted as **g**
_{1} and **g**
_{2}, and Δ**g**̄≡**g**̄** _{1}**-

**g**̄

**. The Hotelling observer computes the linear discriminant**

_{2}where **w**† is the adjoint of the Hotelling template

In practice, numerical computation of **K _{g}**

^{-1}in Eq. (35) can present difficulties due to the large dimensions of

**K**. For example, the 128×128 element detector considered in our simulation studies corresponds to a matrix

_{g}**K**of dimension 16384×16384. To circumvent difficulties with inverting this matrix, a subregion Hotelling observer [37,38] was employed. The subregion Hotelling observer operates on images that have been multiplied by a square window centered at the known signal location, which results in a reduction of the dimensionality of the image autocovariance matrix [37, 38]. In our studies, this was tantamount to extracting a ($\left(\genfrac{}{}{0.1ex}{}{{N}^{2}}{4}\times \genfrac{}{}{0.1ex}{}{{N}^{2}}{4}\right)$) submatrix of

_{g}**K**that was easily inverted and used as a surrogate for

_{g}**K**

_{g}^{-1}in Eq. (35).

A signal-present or signal-absent decision is made by comparing the value of *T*(**g**) to a decision threshold value. When applied to a collection of images, the true positive fraction (TPF) and false positive fraction (FPF) can be estimated (e.g., see Chapter 13 in [23]). In a binary detection task, the TPF and FPF correspond to the probability of detecting the sought-after signal correctly and the false alarm rate, respectively. By varying the decision threshold, a receiver operating characteristic (ROC) curve is obtained [23]. The ROC curve depicts the relationship between the TPF and FPF and completely summarizes the detection performance of the observer.

## 6. Numerical studies

Computer-simulation studies were conducted to investigate the effect of detector spacing on the second-order image statistics and signal detectability in quantitative in-line phase-contrast imaging.

#### 6.1. Computer-simulation studies

A monochromatic and time-harmonic X-ray plane-wave of energy 15 keV was assumed to propagate in the direction of the positive *z*-axis and irradiate an object. A mathematical phantom consisting of uniform ellipsoids that were assigned different values of *n*(*r*⃗) was employed to represent the object. In the detectability study, the object contained only two ellipsoids; a large uniform ellipsoid that represented the object background and a smaller one that represented the signal of interest. The real and imaginary components of the refractive indices were assigned as *δ*=8.6421×10^{-7} and *β*=6.7344×10^{-10} for the background, and *δ*=9.6385×10^{-7} and *β*=1.4236×10^{-9} for the signal. These correspond to values of tumor and fat tissues in human breast. The longest semi-axes for the larger and smaller ellipsoids were 320 *µ*m and 99.2 *µ*m, respectively. Five measurement geometries were considered. In all cases, the first detectorplane was located at *z*
_{1}=9 mm. The second detector-plane (see Fig. 1) was located at *z*
_{2}=37 mm, *z*
_{3}=78 mm, *z*
_{4}=115 mm, *z*
_{5}=176 mm, and *z*
_{6}=250 mm. The detector that contained 128×128 elements of dimension of 5 *µ*m.

From knowledge of the object, the transmitted wavefield *U _{t}* (

*x,y*) on the contact plane was analytically computed according to Eqs. (1)–(3). The 2D fast Fourier transform (FFT) algorithm was employed to compute the 2D convolution of sampled values of

*U*(

_{t}*x,y*) with the Fresnel propagator kernel, and the simulated intensity data

*I*[

_{m}*r, s*] were computed according to Eq. (4). An ensemble containing 40,000 noisy versions of the intensity data was created by generating realizations of uncorrelated Gaussian random process that was characterized by its mean

*µ*and variance

*σ*

^{2}. When generating the noisy data,

*µ*[

*r, s*] was set equal to the noiseless values of

*I*[

_{m}*r, s*] and

*σ*

^{2}=10%. Note that for very low-contrast objects, as considered here, Poisson image noise can also be approximately described as additive and Gaussian.

From the pairs of noisy intensity data corresponding to the different imaging geometries, noisy estimates of *ϕ*̃ * _{m,n}*[

*p,q*] and

*A*̃

*[*

_{m,n}*p,q*] were computed by use of discretized versions of Eqs. (8) and (9) with the 2D FT replaced by the 2D DFT. From these Fourier data, noisy estimates of

*ϕ*[

*r, s*] and

*A*[

*r, s*] were computed by use of the 2D inverse IFFT algorithm. Subsequently, empirical estimates [39] of Cov{

*ϕ*

*[*

_{m,n}*r,s*],

*ϕ*[

_{m,n}*r*′,

*s*′]} and Cov{

*A*[

_{m,n}*r, s*],

*A*[

_{m,n}*r*′,

*s*′]} were computed and compared against the theoretical predictions given by Eq. (31) or (33).

In the signal detection studies, 20,000 noisy estimates of *ϕ _{m,n}*[

*r, s*] and

*A*[

_{m,n}*r, s*] were computed as described above. Half of the images corresponded to signal present and half to signal absent. The subregion Hotelling observer described in Section 5.1 was separately applied to the collections of

*ϕ*[

_{m,n}*r, s*] and

*A*[

_{m,n}*r, s*] images and the TPF and FPF values were computed at 200 different decision threshold values. From the (TPF,FPF) pairs ROC curves were plotted for each imaging geometry.

#### 6.2. Numerical results

**Examples of reconstructed images: **Reconstructed estimates of *ϕ _{m,n}*[

*r, s*] for the different imaging geometries are shown in subfigures (a)–(d) of Fig. 2. From left-to-right, the images correspond to the increasing values of Δ

*z*=

*z*

_{n}-

*z*

_{m}in Fig. 1. The corresponding estimates of

*A*[

_{m,n}*r, s*] are displayed in subfigures (e)–(h). As predicted by the reconstruction formula in Eqs. (8) and (10), the estimates of

*ϕ*[

_{m,n}*r, s*] are contaminated by low frequency noise, which results in a lumpy noisy appearance. As shown in Figs. 2(a)–2(d), this low-frequency noise is to be mitigated when the detector spacing increases. This was also observed by Paganin,

*et al.*, in Reference [17] for the case of the TIE model. Since the reconstruction formula for

*A*[

_{m,n}*r, s*] does not contain any poles, the image estimates contain uncorrelated noise that is not affected by the imaging geometry for the signal-independent noise model considered. These results quantitatively demonstrate that the noise texture in the reconstructed phase and absorption images is fundamentally different.

**Image statistics:**
Figures 3(a) and 3(b) display the images of theoretical and empirical estimates of Cov{*ϕ*
_{1,2}[*r*, *s*],*ϕ*
_{1,2}[0,0]}. The corresponding quantities for Cov{*A*
_{1,2}[*r, s*],*A*
_{1,2}[0,0]} are contained in subfigures (c) and (d). For both the phase and absorption images, the theoretically predicted and empirically determined autocovariance functions appear nearly identical, corroborating the correctness of our analysis. This is more easily confirmed by examination of the profiles through the centers of the autocovariance functions, which are displayed in Fig. 4. As expected, the singularity in the reconstruction formula in Eq. (8) results in strong image correlations in the estimates of *ϕ _{m,n}*[

*r, s*], as shown by Fig. 3(a) or 3(b). Alternatively, Fig. 3(c) or 3(d) demonstrates that the noise in the estimates of

*A*[

_{m,n}*r, s*] is statistically uncorrelated.

Profiles through the autocovariance images corresponding to different detector spacings are shown in Fig. 5. The profiles corresponding to detector spacings Δ*z*=28, 69, 106, 167, and 241 mm are denoted by solid, dashed, dashed-dotted, dotted and thin solid curves, respectively. As predicted by the factor

**Signal detection study:** The ROC curves corresponding to the phase and absorption images for different imaging geometries are shown in Fig. 6.
Figure 6(a) indicates an enhancement in signal detectability in the phase images with increasing values of detector spacing Δ*z*. This is consistent with the noisy appearance of the phase images displayed in Figs. 2(a)–2(d), and the fact that the strength of the noise correlations diminishes with increasing Δ*z*. Note that the ROC curves corresponding to Δ*z*≥106mmnearly overlap with the axes of the figure, which indicates nearly perfect detection performance (due to the relatively low noise level in this example). The case of a higher noise level is considered in Fig. 7, which more clearly demonstrates monotonically increasing detection performance with detector spacing. For the absorption images, the ROC curves in Fig. 6(b) are not strongly affected by the imaging geometry due to the absence of poles in the reconstruction formula and signal-independent noise model. The relatively poor ROC curves are explained by the fact that the absorption contrast of the object is very small. However, a slight enhancement in the object detectability is observed due to the reduction of the image variance with the increased values of the detector spacings.

## 7. Summary and conclusions

In medical imaging applications, task-based measures of image quality are commonly employed to optimize the performance of imaging systems, and there are compelling reasons for employing them elsewhere. A truly meaningful assessment of an imaging system must depend on the task for which the system was intended. Task-based measures of image quality are fundamentally distinct from physical measures in that they require knowledge of the higher-order statistical properties of an image, and take into account the skill of the observer who is interpreting the images. To our knowledge, this is one of the first studies to apply these ideas to quantitative X-ray phase-contrast imaging.

We derived analytic expressions that describe the second-order statistics of the reconstructed absorption and phase images in their continuous and discrete forms. This quantified the spatial correlations that are present in the reconstructed phase images due to the zero-frequency pole in the reconstruction formula, which are absent in the absorption image. The strength of these image correlations was shown to be proportional to the detector-plane spacing. These analyses provided a framework in which the relationships between the second-order image statistics and imaging geometry can be comprehensively understood in quantitative in-line phase-contrast imaging.

Subsequently, knowledge of the images’ covariance matrices and concepts from statistical decision theory were applied to demonstrate how signal detectability in a SKE/BKE detection task was influenced by the image noise statistics corresponding to different imaging geometries. In a similar way, signal detection studies could be conducted to compare or optimize different reconstruction methods and imaging geometries within the context of specific detection tasks.

We assumed the contrast transfer function (CTF) formulation of quantitative phase-contrast imaging. However, our analyses can be readily applied to the TIE formulation. In that case, our qualitative observations regarding the effects of imaging geometry on noise texture and signal detectability will remain unchanged.

In this study we considered the object to be fixed and that randomness in the measurement data was introduced by the detection process only. In a more general study, the object could be treated as random [24] as well as the imaging system itself due to the inherently statistical nature of the partially coherent wavefields that are available in practice. We believe the investigation of these effects on signal detectability in phase-contrast imaging will be important and interesting.

## Acknowledgment

The authors thank Drs. Yongyi Yang and Jovan Brankov for useful discussions of signal detection theory, and benefited greatly from the numerous publications on image quality from Professor Harrison Barrett’s group that inspired this study. This research was supported in part by National Science Council under grant No. NSC 97-2221-E-002-001-MY2 and National Science Foundation Awards CBET-0546113 and CBET-0854430.

## References and links

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

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

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

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

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

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

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

**8. **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(10)**, 2845–2852 (1998). URL http://stacks.iop.org/0031-9155/43/2845. [CrossRef]

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

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

**11. **C.-Y. Chou, Y. Huang, D. Shi, and M. A. Anastasio, “Image reconstruction in quantitative X-ray phase-contrast imaging employing multiple measurements,” Opt. Express **15(16)**, 10,002–10,025 (2007). URL http://www.opticsexpress.org/abstract.cfm?URI=oe-15-16-10002.

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

**13. **T. E. Gureyev and K. A. Nugent, “Phase retrieval with the transport-of-intensity equation. II. Orthogonal series solution for nonuniform illumination,” J. Opt. Soc. Am. A **13(8)**, 1670–1682 (1996). URL http://josaa.osa.org/abstract.cfm?URI=josaa-13-8-1670. [CrossRef]

**14. **A. Barty, K. Nugent, A. Roberts, and D. Paganin, “Quantitative phase tomography,” Opt. Commun. **175(4)**, 329–336 (2000). [CrossRef]

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

**16. **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**, 29,132 (1999). [CrossRef]

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

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

**19. **B. D. Arhatari, A. P. Mancuso, A. G. Peele, and K. A. Nugent, “Phase contrast radiography: Image modelling and optimization,” Rev. Sci. Instrum. **75**, 5271–5276 (2004). [CrossRef]

**20. **X. Wu, H. Liu, and A. Yan, “Optimization of X-ray phase-contrast imaging based on in-line holography,” Nucl. Instrum. Meth. B **234**, 563–572 (2005). [CrossRef]

**21. **Y. Nesterets, S. Wilkins, T. Gureyev, A. Pogany, and A. Stevenson, “On the optimization of experimental parameters for x-ray in-line phase-contrast imaging,” Rev. Sci. Instrum.76(9) (2005).

**22. **T. E. Gureyev, Y. I. Nesterets, A. W. Stevenson, P. R. Miller, A. Pogany, and S. W. Wilkins, “Some simple rules for contrast, signal-to-noise and resolution in in-line x-ray phase-contrast imaging,” Opt. Express **16(5)**, 3223–3241 (2008). URL http://www.opticsexpress.org/abstract.cfm?URI=oe-16-5-3223. [CrossRef]

**23. **H. Barrett and K. Myers, *Foundations of Image Science* (Wiley Series in Pure and Applied Optics, 2004).

**24. **H. H. Barrett, “Objective assessment of image quality: effects of quantum noise and object variability,” J. Opt. Soc. Am. A **7**, 1266–1278 (1990). [CrossRef] [PubMed]

**25. **H. H. Barrett, J. Yao, J. P. Rolland, and K. J. Myers, “Model observers for assessment of image quality,” Proc. Natl. Acad. Sci. **90**, 9758–9765 (1993). [CrossRef] [PubMed]

**26. **K. J. Myers, H. H. Barrett, M. C. Borgstrom, D. D. Patton, and G. W. Seeley, “Effect of noise correlation on detectability of disk signals in medical imaging,” J. Opt. Soc. Am. A **2**, 1752–1759 (1985). [CrossRef] [PubMed]

**27. **C. K. Abbey and M. P. Eckstein, “Classification images for simple detection and discrimination tasks in correlated noise,” J. Opt. Soc. Am. A **24**, B110–B124 (2007). [CrossRef]

**28. **T. E. Gureyev, S. Mayo, S. W. Wilkins, D. Paganin, and A. W. Stevenson, “Quantitative In-Line Phase-Contrast Imaging with Multienergy X Rays,” Phys. Rev. Lett. **86**, 5827–5830 (2001). [CrossRef] [PubMed]

**29. **D. Shi and M. A. Anastasio, “Intensity diffraction tomography with fixed detector plane,” Opt. Eng. **46**, 107,003 (2007). [CrossRef]

**30. **T. E. Gureyev, D. M. Paganin, A. W. Stevenson, S. Mayo, and S. Wilkins, “Generalized eikonal of partially coherent beams and its use in quantitative imaging,” Phys. Rev. Lett. **93(6)**, 068,103 (2004).

**31. **A. Papoulis and S. U. Pillai, *Probability, Random Variables, and Stochastic Processes* (McGraw Hill, 2002).

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

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

**34. **A. R. Pineda and H. H. Barrett, “What does DQE say about lesion detectability in digital radiography?” in Proc. of SPIE , vol. **4320**, pp. 561–569 (2001). [CrossRef]

**35. **R. F. Wagner and D. G. Brown, “Unified SNR analysis of medical imaging systems,” Phys. Med. Biol. **30(6)**, 489–518 (1985). URL http://stacks.iop.org/0031-9155/30/489. [CrossRef]

**36. **W. E. Smith and H. H. Barrett, “Hotelling trace criterion as a figure of merit for the optimization of imaging systems,” J. Opt. Soc. Am. A **3(5)**, 717–725 (1986). URL http://josaa.osa.org/abstract.cfm?URI=josaa-3-5-717. [CrossRef]

**37. **M. Eckstein, J. Bartroff, C. Abbey, J. Whiting, and F. Bochud, “Automated computer evaluation and optimization of image compression of x-ray coronary angiograms for signal known exactly detection tasks,” Opt. Express **11(5)**, 460–475 (2003). URL http://www.opticsexpress.org/abstract.cfm?URI=oe-11-5-460. [CrossRef]

**38. **A. H. Baydush, D. M. Catarious, C. K. Abbey, and C. E. Floyd, “Computer aided detection of masses in mammography using subregion Hotelling observers,” Med. Phys. **30(7)**, 1781–1787 (2003). URL http://link.aip.org/link/?MPH/30/1781/1. [CrossRef]

**39. **M. A. Anastasio, M. Kupinski, and X. Pan, “Noise properties of reconstructed images in ultrasound diffraction tomography,” IEEE Trans. Nucl. Sci. **45**, 2216–2223 (1998). [CrossRef]