Low-noise phase imaging by hybrid uniform and structured illumination transport of intensity equation

Yunhui Zhu; Aamod Shanker; Lei Tian; Laura Waller; George Barbastathis

doi:10.1364/OE.22.026696

1. Introduction

Phase information cannot be obtained directly from electromagnetic waves in the frequency range from THz and above; yet this information is often extremely valuable, as it records the strongest evidence of light interaction with materials of interest. Traditional imaging systems record the intensity of light, which maps to the imaginary part of the complex refractive index n; however, phase provides better contrast for many important classes of objects such as biological cells at visible frequencies, low–Z materials (including liquids and soft tissue at x–ray frequencies, etc.) [1]. The recovery of the real part of n is commonly referred to as “Quantitative Phase Retrieval,” even though the term is slightly misleading—it is not the phase itself that is of interest, but rather the optical path length through the specimen [2]. The latter interpretation has the benefit that it extends the problem naturally to partially coherent systems, where phase is not well defined [3]. Numerous techniques have been developed for the phase retrieval problem, including phase–shifting interferometry [4–6], intensity–based iterative [7, 8] or transport [9] methods, or methods based on Ptychography [10–12].

Among different quantitative phase retrieval methods, the transport of intensity equation (TIE) is convenient because it allows for a simple experimental setup with low demands on both the number and quality of the required optical elements [9, 13]. This feature makes TIE–based methods especially appealing for the X–ray regime; moreover, its demonstrated good tolerance to source incoherence [13–15] also relieves the necessity of using a synchrotron, opening up the possibility of eventual use with tabletop X-ray sources [16].

In the TIE method, the phase of the object is extracted from measurements of the intensity differentials along the propagation axis that occur due to propagation. For pure phase objects, these differentials act as sources and sinks in a Poisson equation which can then be inverted to the phase. However, because the Poisson inversion process significantly amplifies low spatial frequency information, the TIE has an intrinsic problem of amplifying low–frequency noise, which results in cloudy phase results [17,18]. Intuitively, this can be understood because small-scale phase features diffract light significantly, whereas large smooth objects bend the light only slightly, resulting in minimal intensity variations along the propagation direction. The problem is especially destructive in the case of smooth phase objects (e.g. cells, tissue) or low Fresnel numbers (small propagation distances). Many techniques have been proposed to tackle the problem of low frequency noise corruption in TIE, the most popular being computational regularization, linear (Tikhonov) [19] or nonlinear. The latter was recently demonstrated also with sparsity priors, specifically sparsity of hard edges in an otherwise smooth object, to fix the low frequencies [20]. These solutions effectively remove low frequencies, so produce nice results only when low frequency phase information does not exist or is not important. To properly recover low–frequency information, one must change the optical system in order to better capture such information. One way of doing this is to use a large defocus distance in the intensity image, which when combined with small defocus images can produce images with good fidelity across all spatial frequencies [17,18,21–23]. However, large propagation distances may be physically impossible or impractical. Another method places a diffuser in the sample plane in order to amplify the diffraction effects on propagation and thereby achieve better low frequency contrast [24]. Diffraction effects can also be enhanced by structuring the amplitude of the illumination with a spatial light modulator [25,26] and recovering the first gradient of phase for multiple directions, similar to methods in Xray grating interferometry [6, 27]. For example, lateral shear interferomery has been used for phase imaging pure low frequency objects whose high frequency diffraction can be ignored [28, 29]. Recently, Pang et al [30] achieved phase imaging using random illumination patterns with 10 intensity measurements under linearization assumptions.

In this paper, we introduce a novel hybrid TIE quantitative phase retrieval method with optimal response for pure phase objects at all spatial frequency bands. Our phase imaging method only requires three intensity measurements, one with uniform illumination and two with orthogonal sinusoidal (structured) illumination, all at the same defocused plane. By subtracting the diffraction signal with the uniform illumination from that with the structured illumination, we reveal a clean modulation induced signal that is solely depend on the phase gradient. This result is rigorously derived from the original TIE equation without any further approximation or linearization.

Our approach to to low-frequency response enhancement is inspired by structured illumination microscopy [31]. Gustafsson’s insight was that modulating the object intensity shifts the high frequency components of the object down to low frequency region, enabling imaging resolution beyond the detector limit. Essentially, therefore, structured illumination microscopy solves the opposite problem of ours, making up for missing high spatial frequencies to achieve super resolution. Here, we propose to use the same principle, but in reverse: namely, to shift the low frequencies higher, thus making up for TIE’s inherent noise corruption problem. In that sense, it is also analogous to the HiLo reconstruction method [32], which has been shown to be effective in removing low frequency de-focused haze in tissue fluorescence microscopy. Even though in our case the origin of the low-frequency noise is different, we observe similar improvements in the hybrid reconstruction. We also show how to combine and optimize the phase retrieval results from the uniform illuminated TIE and the phase gradient signal to optimally mitigate the effects of the noise at all spatial frequency bands.

The remainder of the paper is structured as follows: in Section 2 we derive the forward and inverse operators for TIE with structured illumination. Based on these, in Section 3 we then design an optimal spectral filter that ideally selects and mixes signal from the uniformly and structure-illuminated TIE images, resulting in a hybrid quantitative phase image with minimum noise power spectrum over the entire spectral domain. Numerical simulation and experimental results with visible light are presented in Section 4.

2. Uniformly and structure-illuminated TIE phase imaging

Let ϕ(x, y) denote the phase of a purely transmitting object located at z = 0, and let I(x, y; d) denote the intensity after propagation by a (small) distance d behind the object or its conjugate in a suitable imaging system. When explicit reference to the lateral coordinates can be omitted without risk of confusion, we will abbreviate as I(d).

The TIE method exploits Fresnel propagation to relate the phase ϕ(x, y) with the intensity distribution differential ∂I(x, y; z)/∂z along the propagation axis z. According to the TIE [9],

{k \frac{\partial I (x, y; z)}{\partial z} |}_{z = 0} = - \nabla_{⊥} \cdot (I (x, y; 0) \nabla_{⊥} ϕ (x, y)),

where x and y are coordinates in the plane perpendicular to the propagation z axis, I and ϕ are the intensity and phase profiles at z = z₀, respectively, k is the wave number of the optical field, and ∇_⊥ ≡ (i∂/∂x + j∂/∂y), is the gradient operator on the x − y plane, with i and j being the unit direction vector in the x and y directions, respectively.

The simplest case is that of a pure phase object and uniform illumination, which is well known but we repeat here for convenience and to establish notation. Let I₀ denote the incident intensity, as shown in Fig. 1(a). The intensity distribution is assumed unchanged (still I₀) right after the thin-phase object. For a pure phase object, Eq. (1) yields

- \frac{1}{k} I_{0} \nabla_{⊥}^{2} ϕ = \frac{\partial I}{\partial z} ≃ \frac{I (d) - I_{0}}{d} \equiv \frac{Δ (d)}{d},

where Δ(d) is the intensity differential signal. This is a Poisson equation and it can be solved in the Fourier domain as

ϕ (x, y) = ℱ^{- 1} [\frac{1}{H (u, v)} ℱ (\frac{Δ (d)}{I_{0}})],

where ℱ is the Fourier transform operator, u and v are the spatial frequency coordinates, and

H (u, v) = \frac{4 π^{2} d}{k} (u^{2} + v^{2})

is the forward transfer function of the Poisson operator. From this result, it is evident that the direct reconstruction described by Eq. (3) is ill-posed at low frequencies and singular at DC—the latter not surprising, since a constant phase factor is usually physically meaningless. Typically, a small regularization parameter is used to avoid the division-by-zero instability, but the spectrum of the intensity differential signal Δ(d) and the noise it carries are still amplified by the inversion process, with an inverse quadratic dependence on spatial frequency. As a result, Eq. (3) typically produces phase estimates that are heavily corrupted by low frequency noise. Figure 2 illustrates this point by plotting the transfer function H and the inverse transfer function H⁻¹ as a function of frequency u.

Fig. 1 Conceptual setup of uniformly (a) and structure-illuminated (b) TIE.

Download Full Size | PDF

Fig. 2 Transfer function (a) and inverse transfer function (b) of uniformly–illuminated TIE (solid line) and structure–illuminated TIE (dashed line). The singular point at spatial frequency u = 0 and low values near low spatial frequency results in noise amplification in this region. The figure also shows increased response and suppressed low–frequency noise amplification at u < f_c (see Eq. (18)), when comparing the structure-illuminated TIE to the uniformly–illuminated TIE.

Download Full Size | PDF

Here, we show that these low–frequency noise problems can be significantly improved by our hybrid structured-illumination TIE method. The idea is that low–frequencies can be shifted to high frequencies through the use of structured illumination. The procedure is shown schematically in Fig. 3. The pure phase object ϕ(x, y) is illuminated with a patterned intensity distribution I_m(x, y, d) as shown in Fig. 1(b). In the simplest case that we investigate here, the modulation is sinusoidal along an arbitrary lateral direction, e.g. x, as

I_{m, X} (x, y, 0) = a I_{0} (1 + η sin (2 π f_{m} x)),

where a is the intensity gain/attenuation, η is the modulation depth, and f_m is a suitably selected spatial frequency. (More on this later.)

Fig. 3 Procedures of phase differential measurement using structure–illuminated TIE.

Download Full Size | PDF

After Fresnel propagation in the modulated illumination case, the intensity differential becomes

\begin{array}{l} k {\frac{\partial I_{m} (x, y; z)}{\partial z} |}_{z = 0} ≃ - k \frac{I_{m} (d) - I_{m} (0)}{d} \equiv - k \frac{Δ_{m} (d)}{d} = \\ = \nabla \cdot (I_{m} (0) \nabla ϕ (x, y)) = \nabla I_{m} \cdot \nabla ϕ (x, y) + I_{m} \nabla^{2} ϕ (x, y), \end{array}

where I_m(d) is the intensity measurement after propagation by distance d. Compared to Eq. (2), an additional cross term ∇I_m · ∇ϕ contributes to the new intensity differential signal Δ_m(d). We can isolate the cross term by combining Eqs. (2) and (6), yielding

\nabla I_{m} \cdot \nabla ϕ = - \frac{k}{d} (Δ_{m} (d) - \frac{I_{m}}{I_{0}} Δ (d)) \equiv - \frac{k}{d} s_{m} (d) .

We refer to s_m(d) as the “modulation cross–term” signal. Note that no movement of the camera is necessary: I₀ and I_m(0) are known (or can be calibrated with a measurement absent any object) while the two measurements I(d) and I_m(d) are obtained with uniform and structured illumination, respectively, downstream at the (out–of–focus) location d.

A derivation similar to Eqs. (5)–(7) has been reported for TIE under general coded illumination [25]; here, we proceed to exploit the equivalence of this result to a convolution (denoted as ⊗) in the spatial frequency domain. Taking the Fourier transform on both sides of this equation, we obtain

ℱ (\nabla I_{m}) \otimes ℱ (\nabla ϕ (x, y)) = - \frac{k}{d} ℱ (s_{m} (d)) .

By the convolution operation, the low-frequency components in the phase differential are shifted into the high–frequency region, where they are immune from TIE’s harmful inverse quadratic amplification. As we noted in Sec. 1, this process is similar to super-resolution imaging using structured illumination microscopy, but in the reverse direction.

Equation (8) is easily solved when I_m is a sinusoidal pattern, as in Eq. (5). The result is

\begin{array}{l} ℱ (ϕ_{x}) (u, v) \otimes π f_{m} a η I_{0} i (δ (u - f_{m}) - δ (u + f_{m})) = \\ = π f_{m} a η I_{0} i (ℱ (ϕ_{x}) (u + f_{m}, v) - ℱ (ϕ_{x}) (u - f_{m}, v)) = - \frac{k}{d} ℱ (s_{m, X} (d)) (u, v), \end{array}

where s_m,X (d) is the modulation cross–term signal after propagation by distance d obtained by using the sinusoidal modulation pattern I_m,X (0). Note that the signal is proportional to the modulation frequency f_m and depth η. The spectrum of the phase differential signal ϕ_x is shifted back to the origin by convolution with the Dirac delta function δ(u + f_m) as

ℱ (ϕ_{x}) (u, v) - ℱ (ϕ_{x}) (u + 2 f_{m}, v) = i \frac{k}{d} \frac{1}{π f_{m} a η I_{0}} ℱ (s_{m, X} (d)) (u + f_{m}, v) .

We can eliminate the double-modulated term ℱ(ϕ_x)(u + 2 f_m, v) in two ways. A simple approach is to use a low pass filter $T_{x} (u, v) = exp [- {(u + f_{m})}^{2} / f_{m}^{2}]$ , which is centered at u = 0 and cut off around u = f_m. This approach is valid if the spectrum of the phase differential has negligible components beyond the modulation frequency f_m, in which case the two lobes of shifted spectrum do not overlap. The phase differential can then be obtained by

ϕ_{x} = ℱ^{- 1} [i \frac{k}{d} \frac{1}{π f_{m} a η I_{0}} ℱ (s_{m, X} (d)) (u + f_{m}, v) \cdot T_{x} (u, v)] .

In principle, the phase profile can be obtained by one more integration as

ϕ = ℱ^{- 1} [H_{m}^{- 1} (u, v) ℱ (\frac{s_{m, X} (d)}{I_{0}}) (u + f_{m}, v) \cdot T_{x} (u, v)],

where

H_{m} (u, v) = \frac{2 π^{2} a η f_{m} d}{k} u

is the forward transfer function for the structured-illumination case. This approach has been used for converting gradient information into phase profile in the grating based X-ray phase contrast imaging techniques [27]. However, the phase retrieval results contain strip-like noise (see below for alternative method and next section for noise analysis). We mention it to emphasize that the transfer function now is linear in the spatial frequency variable where the modulation took place, whereas in the uniform–illuminated case it was quadratic in both spatial frequency variables (compare Eq. (4) and Fig. 2). This explains why the noise influence on low spatial frequencies is expected to be milder in the structured–illumination case, as demonstrated in Fig. 2. The figure also shows that the quadratic response curve eventually catch up and surpass the linear response in the high frequency region beyond some critical frequency f_c.

When the object contains substantial high frequency components beyond f_m, additional measurements can be used to eliminate the double-modulated term ℱ(ϕ_x)(u + 2f_m) [33]. For example, a second modulation cross–term measurement can be obtained by shifting the illumination by π/2

I_{q, X} (x, y, 0) = a I_{0} (1 + η cos (2 π f_{m} x)) .

(Here, “q” stands for “quadrature.”) The resulting modulation cross–term signal s_q,X (d) relates to the spectrum of the phase differential signal as

ℱ (ϕ_{x}) (u, v) + ℱ (ϕ_{x}) (u + 2 f_{m}, v) = - \frac{k}{π f_{m} d a I_{0} η} ℱ (s_{q, X} (d)) (u + f_{m}, v)

Combining Eqs. (10) and (15) we eliminate the double–modulated term and obtain the clean spectrum of the phase differential signal. After inverse Fourier transformation,

ϕ_{x} = ℱ^{- 1} [\frac{k}{d} \frac{1}{2 π f_{m} a η I_{0}} (ℱ (i s_{m, X} (d) - s_{q, X} (d)) (u + f_{m}, v))] .

As before, the phase profile can be obtained in principle by integration; the forward transfer function in this case is

H_{m} (u, v) = \frac{4 π^{2} a η f_{m} d}{k} u

The additional factor of 2 is because we are taking two intensity measurements in this case. Comparing Eqs. (4) and (17), we find that the critical frequency is expressed by

f_{c} = f_{m} / a η .

The direct inversion described in principle in the two above cases (low–pass filtering and quadrature) is problematic, as we will also see in section 4. In practice, we include a second set of measurements using modulation in the perpendicular direction y. Thus, for the quadrature case, we obtain independently the phase differential in the y direction as

ϕ_{y} = ℱ^{- 1} [\frac{k}{d} \frac{1}{2 π f_{m} a η I_{0}} (ℱ (i s_{m, Y} (d) - s_{q, Y} (d)) (u, v + f_{m}))],

where s_m,Y (d) and s_q,Y (d) are modulation cross–term signals with intensity modulation patterns

I_{m, Y} (x, y, 0) = a I_{0} (1 + η sin (2 π f_{m} y))

and

I_{q, Y} (x, y, 0) = a I_{0} (1 + η cos (2 π f_{m} y))

respectively. (The low–pass filter approach can be similarly derived and we omit here.)

Having both ϕ_x and ϕ_y independently available, the most robust way [34, 35] to recover ϕ is by re–constructing the Poisson equation

\nabla^{2} ϕ = \frac{\partial ϕ_{x}}{\partial x} + \frac{\partial ϕ_{y}}{\partial y},

and inverting (similar to Eq. (3)) as

ϕ = ℱ^{- 1} [\frac{1}{4 π^{2} (u^{2} + v^{2})} ℱ (\frac{\partial ϕ_{x}}{\partial x} + \frac{\partial ϕ_{y}}{\partial y})] .

Note that, because Eq. (22) involves taking a derivative, the final solution (23) still maintains inverse linear dependence on the two spatial frequencies u, v. This is beneficial, but further improvement may be obtained by combining the estimate of ϕ obtained from Eq. (23) with the standard Poisson solution (3). This is done in the next section.

3. Error statistics and optimized hybrid TIE reconstruction

Returning to Fig. 2, as we pointed out in the previous section, the inverse transfer function behaves better in the linear case (structured illumination) for spatial frequencies lower than the critical frequency f_c (Eq. (18)) whereas above f_c the reconstruction from uniform illumination should be better. To quantitively describe the “goodness” of imaging strategies, we introduce the power spectral density (PSD) error as a statistical measurement of phase noise corruption. Based on this scheme, we show how to optimally combine the two reconstructions under the assumption of additive white Gaussian noise. The result is our final hybrid TIE reconstruction filter.

To carry out the derivation, we begin by evaluating the noise statistics for four reconstruction strategies, namely:

uniform illumination alone (standard TIE), which we will denote as ϕ̂₁;
sinusoidal illumination along one direction with quadrature recovery, denoted as ϕ̂₂;
sinusoidal illumination alone with two orthogonal direction measurements and quadrature recovery, denoted as ϕ̂₃;
combined uniform and structured illumination, the latter being sinusoidal with two orthogonal direction measurements and quadrature recovery; the exact optimal way of doing the combination will emerge at the end of this section from the preceding analysis, and we will denote the result as ϕ̂_h.

The first three are interesting in their own right, and ϕ̂₁, ϕ̂₃ will be used in the optimal filter derivation. Other reconstruction strategies (e.g., low-pass filter recovery) can also be carried out similarly, but we have omitted them for brevity and because they are sub–optimal in terms of noise performance.

Let the noise–corrupted measurement of the differential signal Δ̂(d) be denoted as

\hat{Δ} (d) \equiv Δ (d) + n (x, y),

where n(x, y) is a wide–sense stationary random process of zero mean. This measurement noise originates from the optical intensity measurement, normally consisting of dark noise of the detector and shot noise from photon statistics. The detector dark noise is uncorrelated for each CCD pixel, and is not dependent on the measured intensity. As a result, the dark noise has a flat spatial power spectrum within the bandwidth of the CCD. The shot noise power scales with the intensity and its spatial spectrum is related to the input intensity profile. However, in typical TIE experiments, the intensity differential Δ(d) and Δ_m(d) are small compared to the initial intensity I₀. For uniform and weakly modulated input fields, the shot noise level of each pixel is also similar to others, resulting in spatially white noise. Under all the above assumptions, the power spectral density of the intensity differential noise n is constant and we define it as

S_{n} (u, v) = σ^{2},

where σ² is the aggregated noise variance.

Starting with the uniformly illuminated case, let ϕ̂₁ denote the noisy estimate of the phase that we obtain from the noisy measurement Eq. (24). Using Eq. (3) and the Wiener–Khintchine theorem, it is straightforward to show that the power–spectral densities S_ϕ(u, v), S_ϕ̂₁ (u, v) are related as

S_{{\hat{ϕ}}_{1}} (u, v) - S_{ϕ} (u, v) \equiv E_{{\hat{ϕ}}_{1}} (u, v) = {| H (u, v) |}^{2} S_{n} (u, v),

where S_n(u, v) is the PSD of n(x, y), and E(u, v) is the PSD error we wish to minimize. Substituting the noise PSD from Eq. (25), the well–known PSD error of the phase retrieval is calculated for uniform TIE in our notation as

E_{{\hat{ϕ}}_{1}} (u, v) = {[\frac{k}{d} \frac{1}{4 π^{2} (u^{2} + v^{2}) I_{0}}]}^{2} σ^{2},

Similarly, we can calculate the noisy phase reconstruction and extract the PSD error of the phase for the sinusoidal-patterned-illuminated TIE. In this case, the PSD error of the phase depends on both the way one recovers phase differentials from Eq. (10) (low–pass filtering or quadrature) and the approach one uses for the final retrieval (direct integral or Poisson reconstruction).

We first consider the PSD error E_{ϕ̂_x} of the phase differential ϕ_x itself. When ϕ_x is obtained by the quadrature method as in Eq. (16), the PSD error becomes

E_{{\hat{ϕ}}_{x}} (u, v) = {(\frac{k}{d} \frac{1}{4 π f_{m} a η I_{0}})}^{2} 4 σ^{2} = {(\frac{k}{d} \frac{1}{2 π f_{m} a η I_{0}})}^{2} σ^{2},

where we also assumed that the noise statistics in the two measurements Δ̂(d) and Δ̂_m(d) are identical. The low–pass filter recovery case yields a similar result, but larger by a factor of 2 (compare Eqs. (13) and (17)).

We next consider the PSD error when relying on a single sinusoidal structured illumination along the x direction, as in Eq. (12), and we denote the reconstruction result as ϕ̂₂. Then

E_{{\hat{ϕ}}_{2}} (u, v) = {(\frac{1}{2 π u})}^{2} E_{{\hat{ϕ}}_{x}} .

As expected, using a second orthogonal sinusoidal illumination, as in Eq. (23), improves the PSD error on the reconstruction ϕ̂₃ as

E_{{\hat{ϕ}}_{3}} (u, v) = {[\frac{1}{4 π^{2} (u^{2} + v^{2})}]}^{2} [{(2 π u)}^{2} E_{{\hat{ϕ}}_{x}} + {(2 π v)}^{2} E_{{\hat{ϕ}}_{y}}] = \frac{1}{4 π^{2} (u^{2} + v^{2})} E_{{\hat{ϕ}}_{x}} .

Again, we have assumed that the PSD error of the phase differential ϕ_y has the same statistics as that of ϕ_x, as is true for systems with orientation symmetry.

Inserting the expression of PSD error of phase differentials from Eq. (28) into the equations above, we finally obtain the corresponding PSD errors of two approaches as

E_{{\hat{ϕ}}_{2}} (u, v) = \frac{1}{u^{2}} {(\frac{k}{d} \frac{1}{4 π^{2} f_{m} a η I_{0}})}^{2} σ^{2}

E_{{\hat{ϕ}}_{3}} (u, v) = \frac{1}{u^{2} + v^{2}} {(\frac{k}{d} \frac{1}{4 π^{2} f_{m} a η I_{0}})}^{2} σ^{2} .

The power spectral densities E_ϕ̂₁ (uniform-illuminated TIE), E_ϕ̂₂ (structured illuminated TIE, forward integral) and E_ϕ̂₂ (structured illuminated TIE, Poisson reconstruction) are shown in Fig. 4 in log scale. We can see, firstly, the clear improvement of noise suppression using two orthogonal sinusoidal illumination steps as in Eq. (23) compared to one-directional integration. This is because the one-directional integration yields an anisotropic noise spectrum—inversely quadratic along the u direction only, and constant along the v direction. As a result, ϕ̂₂ is heavily corrupted by noise on the v axis. The isotropic power spectral density error of E_ϕ̂₃, on the other hand, is inversely linear in all directions. In fact, ϕ₂ is the least-square phase retrieval solution given phase differential profiles in both directions. These conclusions are further demonstrated by the numerical simulations further along in Fig. 5.

Fig. 4 Power spectral density errors of phase retrieval via uniform illuminated TIE (a), sinusoidal-illuminated TIE along one direction (b), sinusoidal-illuminated TIE with two orthogonal directions (c),and the optimized hybrid TIE (d) in logarithmic scale.

Download Full Size | PDF

Fig. 5 Simulation results of phase retrieval using various TIE methods. (a) original phase profile, (b)–(e) phase recovered using uniformly-illuminated TIE ϕ̂₁, structure-illuminated TIE via using direct integral ϕ̂₂ and via Poisson equation construction ϕ̂₃, and hybrid TIE reconstruction ϕ̂_h, (h)–(i) corresponding calculated phase PSD error, plotted in log scale.

Download Full Size | PDF

On the other hand, the isotropic profiles of PSD error spectra E_ϕ̂₁ and E_ϕ̂₃ are consistent with the inverse transfer function behavior depicted in Fig. 2. As predicted in the previous section, noise is relatively suppressed in the low frequency region for the structured illumination TIE but decays relatedly faster in the high frequency region for the uniform illuminated TIE. To obtain the best reconstruction independent of spatial frequency, an optimized phase estimate can be constructed by weighting each of the two signals in the hybrid set of measurements according to its frequency response in different spectral regions. In other words, we can combine high frequency components from the uniformly-illuminated TIE and the low frequency components from the structured-illumination TIE to minimize the overall effect of the noise.

This is formally done as follows: first, we form a convex linear combination of the two reconstructions ϕ̂₁ and ϕ̂₃, filtered respectively as

{\hat{ϕ}}_{h} (x, y) = {\hat{ϕ}}_{1} (x, y) \otimes t (x, y) + {\hat{ϕ}}_{3} (x, y) \otimes (1 - t (x, y)),

where t(x, y) is the impulse response of the shift–invariant filter that we wish to optimize. We assume that t is square integrable and denote its Fourier transform (transfer function) as T(u, v). Forming the autocorrelation of the random process ϕ̂_h and taking its Fourier transform, it is straightforward to show that the PSD error for the new estimate ϕ̂_h is

S_{{\hat{ϕ}}_{h}} (u, v) - S_{ϕ} (u, v) \equiv E_{{\hat{ϕ}}_{h}} (u, v) = {| T (u, v) |}^{2} E_{{\hat{ϕ}}_{1}} + {| 1 - T (u, v) |}^{2} E_{{\hat{ϕ}}_{3}} .

Now we assume the signals are band limited, and exclude from any consideration a small region of spatial frequencies around (u, v) = (0, 0). The latter exclusion is justified because the DC component of a phase signal has no physical meaning and so we are justified in ignoring it. Under these assumptions, the total noise power P_N can be computed as

P_{n} = \iint_{★} [{| T (u, v) |}^{2} E_{{\hat{ϕ}}_{1}} + {| 1 - T (u, v) |}^{2} E_{{\hat{ϕ}}_{3}}] d u d v,

where ∬_★ denotes the integral over all spatial frequencies excluding an appropriately small region around the origin, as noted just before.

To minimize P_n, we use the variational method. By setting the variation of P_n to equal zero,

δ P_{n} = \iint 2 (T E_{{\hat{ϕ}}_{1}} + (T - 1) E_{{\hat{ϕ}}_{3}}) δ T d u d v = 0,

and substituting Eqs. (27) and (30), we obtain the optimal filter transfer function

T (u, v) = \frac{E_{{\hat{ϕ}}_{3}} (u, v) (u, v)}{E_{{\hat{ϕ}}_{1}} + E_{{\hat{ϕ}}_{3}} (u, v)} = \frac{u^{2} + v^{2}}{{(f_{m} a η)}^{2} + (u^{2} + v^{2})} .

Interestingly, this filter is a Lorentzian function with a half width equal to the critical frequency f_c, and we refer to it as the hybrid TIE (hybrid uniformly and structured illuminated TIE) filter. It yields the optimal (in the sense of minimal overall noise power) PSD error

E_{{\hat{ϕ}}_{h}} (u, v) = \frac{E_{{\hat{ϕ}}_{1}} (u, v) E_{{\hat{ϕ}}_{3}} (u, v)}{E_{{\hat{ϕ}}_{1}} (u, v) + E_{{\hat{ϕ}}_{3}} (u, v)} .

The minimized noise spectrum is shown Fig. 4(d). We can see that, indeed, the overall noise spectrum is reduced compared to all the power spectra E_ϕ̂₁, E_ϕ̂₂ and E_ϕ̂₃.

4. Simulation and experimental results

Here we report simulation and experimental results to confirm and characterize the performance of various structured illumination approaches for TIE phase retrieval. We first simulate the phase imaging in the x-ray region. A pure phase profile of a 2×2 mm² Lena’s image, as shown in Fig. 5(a), is used as the object. The largest phase value is set to ±2.6 rad, while the background is 0 rad. We assume that the illumination wavelength is in the hard x-ray range of 0.0633 nm (∼20 keV), propagation distance d = 0.4 m and the sinusoidally patterned illumination has a single spatial frequency f_m at 14.2 mm⁻¹ with 50% modulation depth. Pixel size is assumed to be 4.4 μm and space–bandwidth product is set to 1024×1024. Additive white noise in the intensity differential measurement is assumed to be σ² = 10⁻⁴, normalized to background intensity power $I_{0}^{2} = 1$ . All these numbers are within practical range of realistic x-ray experiment. Note that the propagation distance of 0.4 m is far smaller than the Talbot distance of the grid, thereby is insufficient for obtaining diffraction signal from low frequency components (< f_m) using uniform illumination TIE. However, multiplication of the propagation distance quickly becomes impractical for in-door operations. We show that structure illumination can revive low frequency signals under propagation distance constrain.

Figure 5(b) shows the phase reconstruction using uniformly-illuminated TIE, corrupted with low-frequency noise. The structured illuminated TIE phase retrieval ϕ̂₂ shown in Fig. 5(c) contains substantial strip-like high–frequency noise along the x direction, as predicted (see Fig. 4(b).) The low frequency background noise is substantially suppressed for ϕ̂₃ as shown in Fig. 5(d), where we used the modulation in both orthogonal directions. However, the image is not as sharp as in ϕ̂₁. The hybrid phase reconstruction ϕ̂_h shown in Fig. 5(e) both removes noise effectively and preserves sharp and smooth features alike, resulting in a high fidelity phase retrieval. In Figs. 5(f)–5(i) we also show the corresponding PSD error, as computed from the reconstructions, showing good agreement with the theoretical predictions given in Fig. 4.

We next built a TIE phase imaging experiment in the optical range. The experiment is designed to mimic the x-ray simulation presented above. We use the operation wavelength of 633 nm and a propagation distance of 0.4 mm, while keeping the pixel size at 4.4 μm. By keeping the Fresnel number (dx)²/λd fixed, the TIE signal stays the same. We use a lens (N-BK7 Half-Ball Lens 6.0 mm Diameter from Edmund Optics) as the phase object. To reduce the total phase difference thereby improving the condition of the inverse problem, the lens object was immersed in index-matching oil. At the operation wavelength of 633 nm, the refractive index difference between the lens and the immersion oil is effectively reduced to 0.0036±0.0001, yielding a nominal maximum phase variation of 106.6 rad. We chose the lens object deliberately because (i) its shape is calibrated so it provides ground truth and (ii) it contains mostly low spatial frequency features, thus it accentuates TIE’s problems at that band.

The experimental arrangement is shown in Fig. 6. A multimode fiber delivers the light beam from a broadband LED illuminator to a focusing lens. The beam is focused down to a pinhole (diameter 50 μm) and then collimated to provide quasi–plane wave illumination. A narrow band (1 nm @ 633 nm) spectral filter is inserted in the collimator. In this way, we obtain uniform illumination with high spatial and temporal coherence. The object is imaged with a magnification ×0.5 via a 4f telescope system onto the detector. The lenses used in the telescope system are 2-inch aspheric lenses with focal length of 20 cm and 10 cm, respectively. The detector has pixel size of 4.4 μm. It is placed on a rail that can be translated along the optical propagation axis z. In the uniformly-illuminated TIE experiment, we move the detector to obtain images both in focus and at a defocus distance d = 0.1 mm, corresponding to a propagation distance of d = 0.4 mm for the object. The additional in-focal imaging is taken for background normalization. For the structured illuminated TIE, we only take the defocused images with the patterned illumination, which is generated by inserting a sinusoidal amplitude grating (spatial frequency 10 mm⁻¹, modulation depth of 40%, Applied Imaging Inc.) immediately before the lens, but within the depth of focus such that it is effectively collocated with the sample. Both horizontal and vertical structured illumination were implemented in the experiment and one defocused image is taken for each. The phase differential is obtained via the filter in Eq. (11) and the phase profile is retrieved via the Poisson Eq. (23). Note that the 4-f system is only used to separate the object and image plane while keeping the effective propagation distance to 0.4 mm. No operation on the optical Fourier plane is implemented, as no 4-f system is available in the x-ray region.

Fig. 6 Experimental setup for structure-illuminated TIE. The object (a lens) is placed in approximately the same plane as a sinusoidal amplitude grating for obtaining the illumination structure.

Download Full Size | PDF

Figure 7 shows the experimental phase retrieval results. As expected due to the low frequency content, the uniformly-illuminated TIE result (Fig. 7(a)) is gravely corrupted by noise. On the other hand, the structured–illumination approach on the other hand, successfully captures the phase profile (Fig. 7(b)), but misses the sharp boundary. Finally, the hybrid reconstruction ϕ_h shown in Fig. 7(c), also obtains the low–frequency phase profile successfully. Moreover, from Fig. 7(c) and the cross–sections shown in Fig. 7(d), we observe that a small bubble in the immersion oil was accidentally trapped in the vicinity of the lens. This bubble is not visible in Fig. 7(c) where, again as expected, high frequencies/small features are filtered out. The overall measured SNR (bubbles ignored) for the uniform, structure illumination TIE and the hybrid reconstruction are 0.79, 4.51 and 4.77, respectively. It is noted that the noise spectrum for the uniform TIE result shown in Fig. 7(a) is heavily concentrated in the low frequency domain. This could be the result of intensity noise due to propagation of illumination non-uniformities. Further optimization of the hybrid reconstruction can be made based on the noise spectrum obtained in the experiment.

Fig. 7 Phase retrieval of (a) (failed) uniformly-, (b) structure-illuminated TIE, and (c) the hybrid reconstruction of a lens. (d) shows an extracted line at y = 2.112 mm. Note that an air bubble is visible in the direct TIE and hybrid reconstruction, but not captured in the structured illuminated case. The bubble area is labeled with a small black box in (b) and (c).

Download Full Size | PDF

5. Conclusion

We have demonstrated both theoretically and experimentally that substantial noise suppression is obtained, especially in the low–frequency region, by a hybrid reconstruction of TIE phase retrieval using both uniform and structured illumination. Future work will focus on designing the illumination pattern optimally for a given situation [25], as well as exploring the possibility of also achieving enhanced resolution through proper design on of the pattern. The latter has been done recently for phase imaging with interferometry [36] and can be incorporated into our algorithm in the future. It is noted that so far we have not utilized sparsity or other compressive priors. These should be expected to further suppress noise and improve reconstruction quality.

This research was funded by the US Department of Homeland Security. We would like to thank Zhengyun Zhang for helpful discussions.

References and links

1. P. Bach, J. Jett, U. Pastorino, M. Tockman, S. Swensen, and C. Begg, “Computed tomography screening and lung cancer outcomes,” Jama 297(9), 953–961 (2007). [CrossRef] [PubMed]

2. K. A. Nugent, D. Paganin, and T. E. Gureyev, “A phase odyssey,” Phys. Today 54(8), 27–32 (2001). [CrossRef]

3. J. W. Goodman, Statistical Optics (Wiley-Interscience, 1985).

4. R. K. Crane, “Interference phase measurement,” Appl. Opt. 8, 538 (1969).

5. O. Y. Kwon, “Multichannel phase-shifted interferometer,” Opt. Lett. 9(2), 59–61 (1984). [CrossRef] [PubMed]

6. T. Weitkamp, A. Diaz, C. David, F. Pfeiffer, M. Stampanoni, P. Cloetens, and E. Ziegler, “X-ray phase imaging with a grating interferometer,” Opt. Express 13(16), 6296–6304 (2005). [CrossRef] [PubMed]

7. R. Gerchberg and W. Saxton, “Phase determination for image and diffraction plane pictures in the electron microscope,” Optik 34(3), 275–284 (1971).

8. J. R. Fienup, “Reconstruction of an object from the modulus of its Fourier transform,” Opt. Lett. 3(1), 27–29 (1978). [CrossRef] [PubMed]

9. M. R. Teague, “Deterministic phase retrieval: a Green’s function solution,” J. Opt. Soc. Am. 73(11), 1434–1441 (1983). [CrossRef]

10. J. M. Rodenburg and H. M. Faulkner, “A phase retrieval algorithm for shifting illumination,” Appl. Phys. Lett. 85(20), 4795–4797 (2004). [CrossRef]

11. P. Thibault, M. Dierolf, O. Bunk, A. Menzel, and F. Pfeiffer, “Probe retrieval in ptychographic coherent diffractive imaging,” Ultramicroscopy 109(4), 338–343 (2009). [CrossRef] [PubMed]

12. G. Zheng, R. Horstmeyer, and C. Yang, “Wide-field, high-resolution Fourier ptychographic microscopy,” Nat. Photon. 7(9), 739–745 (2013). [CrossRef]

13. D. Paganin and K. A. Nugent, “Noninterferometric phase imaging with partially coherent light,” Phys. Rev. Lett. 80(12), 2586 (1998). [CrossRef]

14. T. E. Gureyev, A. Roberts, and K. A. Nugent, “Partially coherent fields, the transport-of-intensity equation, and phase uniqueness,” J. Opt. Soc. Am. 12(9), 1942–1946 (1995). [CrossRef]

15. 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(2), 569–580 (2006). [CrossRef]

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

17. L. Waller, L. Tian, and G. Barbastathis, “Transport of intensity phase-amplitude imaging with higher order intensity derivatives,” Opt. Express 18(12), 12552–12561 (2010). [CrossRef] [PubMed]

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

19. L. Turner, B. Dhal, J. Hayes, A. Mancuso, K. Nugent, D. Paterson, R. Scholten, C. Tran, and A. Peele, “X-ray phase imaging: Demonstration of extended conditions for homogeneous objects,” Opt. Express 12(13), 2960–2965 (2004). [CrossRef] [PubMed]

20. L. Tian, J. C. Petruccelli, and G. Barbastathis, “Nonlinear diffusion regularization for transport of intensity phase imaging,” Opt. Lett. 37(19), 4131–4133 (2012). [CrossRef] [PubMed]

21. E. D. Barone-Nugent, A. N. T. O. N. Barty, and K. A. Nugent, “Quantitative phase–amplitude microscopy I: optical microscopy,” J. Microsc. 206(3), 194–203 (2002). [CrossRef]

22. C. Zuo, Q. Chen, Y. Yu, and A. Asundi, “Transport-of-intensity phase imaging using Savitzky-Golay differentiation filter-theory and applications,” Opt. Express 21, 5346–5362 (2013). [CrossRef] [PubMed]

23. Z. Jingshan, R. A. Claus, J. Dauwels, L. Tian, and L. Waller, “Transport of Intensity phase imaging by intensity spectrum fitting of exponentially spaced defocus planes,” Opt. Express 22(9), 10661–10674 (2014) [CrossRef] [PubMed]

24. P. F. Almoro, L. Waller, M. Agour, C. Falldorf, G. Pedrini, W. Osten, and S. G. Hanson, “Enhanced deterministic phase retrieval using a partially developed speckle field,” Opt. Lett. 37(11), 2088–2090 (2012). [CrossRef] [PubMed]

25. A. Shanker, L. Tian, and L. Waller, “Defocus-based quantitative phase imaging by coded illumination,” Proc. SPIE 8949, 89490R (2014).

26. Y. Zhu, A. Pan, and G. Barbastathis, “Low-noise TIE phase imaging by structured illumination,” in Computational Optical Sensing and Imaging, Optical Society of America (2014), CTh3C-5. [CrossRef]

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

28. K. Ichikawa, A. W. Lohmann, and M. Takeda, “Phase retrieval based on the irradiance transport equation and the Fourier transform method: experiments,” Appl. Opt. 27, 3433–3436 (1988). [CrossRef] [PubMed]

29. F. Roddier, “Wavefront sensing and the irradiance transport equation,” Appl. Opt. 29, 1402–1403 (1990). [CrossRef] [PubMed]

30. P. Gao, G. Pedrini, C. Zuo, and W. Osten, “Phase retrieval using spatially modulated illumination,” Opt. Lett. 39, 3615–3618 (2014). [CrossRef] [PubMed]

31. M. G. Gustafsson, “Nonlinear structured-illumination microscopy: wide-field fluorescence imaging with theoretically unlimited resolution,” Proc. Natl. Acad. Sci. U.S.A. 102(37), 13081–13086 (2005). [CrossRef] [PubMed]

32. J. Mertz and J. Kim, “Scanning light-sheet microscopy in the whole mouse brain with HiLo background rejection,” J. Biomed Opt. 15(1), 016027 (2010). [CrossRef] [PubMed]

33. M. R. Arnison, K. G. Larkin, C. J. Sheppard, N. I. Smith, and C. J. Cogswell, “Linear phase imaging using differential interference contrast microscopy,” J. Microsc. 214(1), 7–12 (2004). [CrossRef] [PubMed]

34. S. Velghe, J. Primot, N. Guérineau, M. Cohen, and B. Wattellier, “Wave-front reconstruction from multidirectional phase derivatives generated by multilateral shearing interferometers,” Opt. Lett. 30(3), 245–247 (2005). [CrossRef] [PubMed]

35. S. Ettl, J. Kaminski, M. C. Knauer, and G. Häusler, “Shape reconstruction from gradient data,” Appl. Opt. 47(12), 2091–2097 (2008). [CrossRef] [PubMed]

36. S. Chowdhury and J. Izatt, “Structured illumination quantitative phase microscopy for enhanced resolution amplitude and phase imaging,” Biomed. Opt. Express 4(10), 1795–1805 (2013). [CrossRef] [PubMed]

Low-noise phase imaging by hybrid uniform and structured illumination transport of intensity equation

Abstract

1. Introduction

2. Uniformly and structure-illuminated TIE phase imaging

3. Error statistics and optimized hybrid TIE reconstruction

4. Simulation and experimental results

5. Conclusion

References and links

Cited By

Figures (7)

Equations (38)

Optics Express