1. Introduction

The transport of intensity equation (TIE) is a second order elliptical, non-separable and inhomogeneous partial differential equation which relates the irradiance and the variation of the irradiance along the direction of propagation to a Laplacian-like function of the phase:

The symbol $\overrightarrow{\nabla }$ denotes the gradient in the plane normal to the beam direction, I is the image intensity at the plane of interest normal to the optic axis at z, φ is the phase change induced by the sample, k is the wave vector and z is the coordinate along the direction of propagation. Since the right-hand side, ∂I/∂z, is accessible through experiment as a finite difference derivative in the defocus, the phase shift, φ, can in principle, be retrieved by solving the TIE. Various authors have proposed solutions to this problem [1–4], see [5] for an overview.

For the TIE to have a unique solution, the boundary conditions of the problem need to be specified. These are often taken as periodic [6], or Dirichlet or Neumann boundary conditions are assumed on the measurements’ edges [3] or in a user-defined region [7]. As it is often difficult to impose realistic boundary conditions, one might opt to go without them and instead consider the TIE an underdetermined problem and remedy its ill-conditionedness and non-uniqueness by imposing one or more additional constraints to the solution. A limited degree of spatial coherence often causes the additional problem in many experimental setups that the low spatial frequency information of the $\frac{\partial I}{\partial z}$ measurement is either noisy (e.g. if the defocus has been changed only over a limited range, or at spatial frequencies approaching the lateral coherence length) and/or contains systematic errors (e.g. by some incoherent background signal or the neglect of image distortions when obtaining this measurement by a finite difference approach). Both, problems, unspecified boundary conditions, as well as unreliable low-frequency components of the measurement can be addressed by introducing physically reasonable constraints to the solution.

Oftentimes, the sought-after solution is the one sparsest in a basis suitable for the problem at hand. A naive implementation would then attempt to minimize the ℓ₀-norm of the solution in said basis. However, the non-convexity of the ℓ₀ constraint makes gradient-based optimization nearly impossible. It is the compressed sensing (CS) community’s great achievement to have shown that under general and reasonable circumstances, minimization of the convex ℓ ₁-norm leads to the sparsest solution all the same, thus making the problem treatable with gradient-based optimization techniques.

Total variation (TV) regularization, which aims at minimizing the ℓ₁-norm of the x- and y-derivatives of the solution, (see Eq. 3) seems to have been established as a physically reasonable constraint for reducing low spatial frequency artifacts in TIE phase retrieval [8, 9]. Since the resulting solutions exhibit sparse derivatives, TV-regularization should only be considered for piece-wise constant objects.

In this work a variety of example problems is presented where conventional TV-regularization falls short and where, as expected, especially the lower spatial frequencies are problematic. It is then demonstrated that a substantial improvement in reconstruction quality is obtained if instead the ℓ₀-norm of the x- and y-derivatives, symbolically denoted ∥∇∥₀ in this work, is driven to below a certain user-defined threshold.

Because of the discontinuous nature of the ℓ₀ norm, gradient-based optimization cannot be invoked to make the solutions fulfill this non-convex constraint; hence the gradient flipping algorithm (GFA) is applied. GFA has been inspired by a class of algorithms known in crystallography as charge flipping (CF) algorithms [10], which have been shown to optimize problems with non-convex constraint sets and are able to overcome stagnation at local optima [11].

Since in both TV-regularization and ∥∇∥₀ regularization, the basis which is constrained to be sparse, the gradient ∇ is the same, the improvements presented in this paper may solely be attributed to GFA’s ability to deal with the non-convex constraint imposed by ∥∇∥₀ regularization.

In Section 2 the principle of TVAL3, a state-of-the-art TV-constrained PDE solver [12, 13], are summarized, and GFA is explained. In Section 3.1 the convergence of GFA is verified by comparing reconstructions of piece-wise linear objects to those obtained by TVAL3. In Section 3.3 the reconstruction of an object that’s only partially piece-wise constant, TV-regularization contains significant low spatial frequency artifacts, ∥∇∥₀ yields a much improved result. The results in Section 3.4 show ∥∇∥₀’s superiority over TV-regularization when piece-wise linear objects are to be reconstructed. In Section 3.5 the robustness of ∥∇∥₀-regularization in the presence of physical wave propagation is shown. These assertions bear out for our test on experimental data presented in Section 4. In Section 5 the conclusions are drawn.

2. Reconstruction algorithms

In the following section, we explain the two distinct schemes to constrain the solution of the TIE, namely total variation minimization by the augmented Lagrangian method [14] and alternating direction [15] algorithms (TVAL3) and the gradient flipping algorithm (GFA) [5].

In the remainder of this paper (except for the experimental case presented in section 4), pure phase objects are assumed, leading to in-focus images of uniform intensity. Equation 1 thus simplifies to

2.1. TVAL3

In the TVAL3 scheme a compressed sensing problem with TV regularization is considered. Rewritten for the notations used in this paper, it reads

In TVAL3, this problem is solved by minimizing the augmented Lagrangian, which is a combination of a squared penalty function and the classic Lagrangian function [13, 14]. This alleviates difficulties in connection to large squared penalty terms such as ill-conditioning and bias. As shown in [12], Eq. 3 can be solved fast and efficiently with an alternating directions approach.

It is worth mentioning that C. Li compared the performance of the TVAL3 with other packages such as l₁-magic [16, 17], TwIST [18] and NESTA [19] for different scenarios [13]. It was demonstrated that TVAL3 outperforms other stat-of-art implementations and has the potential to solve compressed sensing problems with TV minimization in an affordable time with high accuracy [12, 13].

Since TVAL3 is open source, it can be adjusted to one’s specific needs. For this paper, the measurement matrix, ∇², was implemented implicitly as a two-dimensional convolution with the kernel

For the reconstruction from experimental data in Sec. 4, absorption needs to be accounted for and thus the operator ∇I∇ from Eq. (1) is implemented implicitly in TVAL3. This was done as a three-stage process, first the x- and y-derivatives are approximated as a central finite difference of second order accuracy, then both derivatives are multiplied with I, and finally both the x-and y-derivatives are derived with respect to x and y, respectively, and summed up. Again, the need for explicit boundary conditions was circumvented by not computing the Laplacian at the edges of the image, so that an image ϕ with dimensions N × N, yields ∇(I∇ϕ) with dimensions (N − 4) × (N − 4). The estimated phase image is now four pixels larger in both directions than the measurements.

2.2. Gradient flipping algorithm

In this section the gradient flipping algorithm (GFA) developed to solve the ∥∇∥₀ problem will be described. The general expression Eq. (1) of the TIE is used to formulate the GFA. Reformulation for the special case Eq. (2) is trivial as one simply has to set I = 1. The numerical experiments in Section 3 all assume the special case Eq. (2) and thus no unfair advantage was given to ∥∇∥₀ over TV-minimization.

The GFA is introduced to solve the ∥∇∥₀-problem of finding the phase φ that has a certain fraction ε of its gradients equal to zero whilst satisfying the TIE, i.e.

The GFA makes use of the following two quantities, obtained from rearranging Eq. (1),

At the ℓ^th iteration of the GFA, first the gradients are calculated

The free parameter R_LP can be chosen to minimize

The iterations start at ℓ = 1 and D⁽⁰⁾ is set to $d{I}_z^{\exp }$. When the algorithm is stopped, at iteration ℓ, the phase φ is calculated from

The gradient flipping step in Eq. (9) is analogous to that in the charge flipping algorithms [10, 11] in crystallography, where the somewhat counter-intuitive notion was used that flipping the sign of entries with values below a certain positive threshold, drives these values to zero. In this work the same notion is employed to reach the goal in Eq. (5), with the small modification that the absolute value of the entries is compared against the threshold. This class of algorithms has proven to optimize problems with non-convex constraint sets and to be able to overcome stagnation at local optima [11].

3. Numerical experiments

Choosing various numerical phase-objects, we aim to investigate the performance of both algorithms namely, TVAL3 as well as GFA, for different scenarios namely, periodic and non-periodic boundary condition, partially piece-wise, piece-wise linear, and the effect of noise and physical wave propagation. For all different scenarios, the intensity of the principal (focused) plane is assumed to be unity over the entire plane to mimic the pure phase-object condition. With the exception of Sec. 3.5, the right hand-side of Eq. (1) is computed by passing the original phase to the implemented two dimensional Laplace operator which resembles the ideal intensity variation along the optical axis.

Since these are numerical experiments, there is access to the object’s phase and hence the root-mean-square-error (RMSE) of the reconstructions can be computed. In order to compare TVAL3 and GFA under their best possible performances, in this section, their respective free parameters (μ and β for TVAL3, δ and R_LP for GFA) are manually adjusted to minimize the RMSE.

3.1. Validation of GFA for piece-wise constant phase objects

A piece-wise constant head-phantom obtained by the Matlab expression ‘phantom()’ is taken to be a pure phase-object with periodic boundary conditions as illustrated in Fig. 1(a). The intensity variation is shown in Fig. 1(d) and the reconstructed phase by means of GFA and TVAL3 are depicted in Figs. 1(b) and 1(c), respectively. The free parameters of the TVAL3 were optimized in such a way that no further improvement of the result could be observed. Moreover, a line profile of the data was extracted along the red line.

 

Fig. 1 Periodic piece-wise constant phase object: a) Head-phantom original phase. b) Reconstructed phase by GFA. c) Reconstructed phase by TVAL3 d) Intensity variation along the optical axis. e) Phase cross sections taken along the red-line in (a).

Download Full Size | PDF

In [7] problems have been pointed out with conventional Fourier based reconstruction methods in the absence of periodic boundary conditions, even when symmetrization was applied. In order to investigate the reconstruction ability of the algorithms in the case of non-periodic phase-objects, a shifted head-phantom depicted in Fig. 2(a) has been considered. The phases reconstructed by means of GFA and TVAL3 are shown in Figs. 2(b) and 2(c), respectively. The Laplacian of the phase is shown in Fig. 2(d). Both algorithms retrieve the phase information faithfully. The slightly lower resolution associated with GFA is addressed in Sec. 3.5, where both methods turn out to have the same resolution equal to that of the measurements when they stem from the physical process of Fresnel propagation.

 

Fig. 2 Non-periodic piece-wise constant phase object: a) Head-phantom original phase. b) Reconstructed phase by GFA. c) Reconstructed phase by TVAL3. d) Laplacian of the phase. e) Phase cross sections taken along the red-line in (a).

Download Full Size | PDF

3.2. Reconstruction from noisy images

In order to evaluate the capability of both algorithms to reconstruct the phase information in the presence of noise, we added 10% and 30% Gaussian noise to the Laplacian of the phases. Figures 3(a) and 3(b) show the reconstruction result of the head-phantom by means of GFA. Furthermore, Figs. 3(c) and 3(d) illustrate the retrieved phase information employing the TVAL3 algorithm. The original phase is depicted in Fig. 3(e).

 

Fig. 3 a) Reconstruction in the presence of 10dB noise by GFA. b) Reconstruction in the presence of 30dB noise by GFA. c) Reconstruction in the presence of 10dB noise by TVAL3. d) Reconstruction in the presence of 30dB noise by TVAL3. e) Original phase.

Download Full Size | PDF

Moreover, also an FFT-based solver using Tikhonov regularization, as described in [20], is employed on this data, since it is a standard method used in the presence of noise. In order to choose the regularization parameter α optimally, following the procedure for GFA and TVAL3, the RMSE was calculated for reconstructions with various α-values and Figs. 4(a) and 4(b) show the phases reconstructed with the optimum value.

 

Fig. 4 Results of the Tikhonov-regularized FFT-solver. a) Reconstructed phase for 10 dB noise, RMSE = 10.6. d) Reconstruction for 30 dB noise, RMSE = 11.4.

Download Full Size | PDF

In order to compare the quality of the reconstructed phases, RMSE is used as a figure of merit. Table 1 shows the RMSE for the aforementioned scenarios.

Table 1. RMSE of the TV-regularization, ∥∇∥₀-regularization and Tikhonov regularization as applied to the various test-cases in this paper. ^*Periodic boundary conditions. ^**Non-periodic boundary conditions.

View Table

3.3. Partially piece-wise constant phase objects

To further investigate the performance of the GFA as well as the TVAL3 algorithm for a proper recovery of low frequency information, we constructed the partially piece-wise constant phase object shown in Fig. 5(a). This phase map can be obtained by the Matlab expression ‘membrane()’. Figures 5(b) and 5(c) illustrates the reconstructed phase by means of GFA as well as TVAL3, respectively. Figure 5(d) depicts the intensity variation $d{I}_z^{\exp }$. It is worth noting that the TVAL3 parameters were optimized in multiple trials to obtain the best possible reconstruction, i.e. the one with the lowest RMSE. Moreover, the threshold parameter δ in GFA is refined during the reconstruction such that the number of flipped pixels is equal to one quarter of the field of view. The RMSE values of the reconstructed phase maps are shown in Table 1.

 

Fig. 5 Partially piece-wise constant phase object: a) Original phase. b) Reconstructed phase by GFA (RMSE = 0.22). c) Reconstructed phase by TVAL3 (RMSE = 5.4) d) Intensity variation in transverse plane.

Download Full Size | PDF

3.4. Piece-wise linear phase objects

Finally, the cameraman phantom shown in Fig. 6(a) from the Matlab built-in library of demo images was chosen as a piece-wise linear object to investigate and evaluate the capability of both phase retrieval algorithms to deal with those types of objects. Figure 6(b) depicts the reconstruction phase by means of the TVAL3 approach in which μ = 2⁵ and β = 2⁴. Figure 6(c) illustrates the retrieved phase information by means of the GFA approach. The measurement matrix obtained by applying (4) to the original phase is shown graphically in Fig. 6(d). It worth noting that again, for the TVAL3 based reconstruction the parameters have been chosen to minimize the RSME (see Table 1).

 

Fig. 6 Piece-wise linear phase object: a) Original phase. b) Reconstructed phase by TVAL3 (RMSE = 1.7). c) Reconstructed phase by GFA (RMSE = 0.95) d) Graphical representation of the measurement.

Download Full Size | PDF

3.5. Fresnel-propagation based measurements

In any practical experiment, the quantity $d{I}_z^{\exp }$ is calculated through a finite difference method from a series of defocused images, and this is only an imperfect approximation of the true Laplacian. In this section the influence of this approximation is tested by generating two images of the head-phantom phase object through Fresnel-propagation at over- and under-focus and taking $d{I}_z^{\exp }$ as the finite difference of both images.

By setting the wavelength of the illumination to 0.5 μm, the pixel size to 1 μm and the numerical aperture to 0.125, it was ensured that the images are sampled at half the Nyquist frequency, as is often the case in practice. The over- and under-focus values were set to 0.1 μm and −0.1 μm respectively to approximate the Laplacian as good as possible; it is noted that in practice much larger values are used (for instance, see Sec. 4) which generally leads to a worse approximation. The Fresnel-propagation based measurements are depicted in Fig. 7(a).

 

Fig. 7 Fresnel-propagation based measurement. a) Finite difference of two images at 0.1 μm over- and under-focus. b) Phase reconstructed by GFA (RMSE = 1.50). c) Phase Reconstructed by TVAL3 (RMSE = 1.95). d) Plot of the radially averaged power spectrum. e) Line profile of the data extracted along the red line.

Download Full Size | PDF

Figures 7(b) and 7(c) depict the phases reconstructed by means of GFA and TVAL3, respectively. The low frequency artifacts are quite evident in the latter and are reflected in the respective RMSEs of 1.50 and 1.95, see Tab. 1. In order to present a more detailed comparison, Fig. 7(d) shows the radially averaged power spectrum of the measurement and the reconstructed phases, clearly demonstrating that the resolution of the retrieved phases is predominantly determined by that of the measurements, i.e. the numerical aperture. Furthermore, Fig. 7(e) shows the line profile of the data extracted along the red line.

4. Experiment

Contrary to Sec. 3, no access to the real phases is available. The free parameters of the reconstruction algorithms (μ and β for TVAL3, δ and R_LP for GFA) are adjusted manually to obtain a minimum χ², as defined in Eq. (13). For Tikhonov regularization FFT-based the α-value from [20] has been adjusted manually to yield a realistic dynamic range in the reconstruction without displaying the low frequency artifacts typical for non-regularized FFT solutions.

Illustrated graphically in Fig. 8, a simple optical setup is employed to evaluate the performance of the two algorithms. The setup is comprised of a green laser emitting coherent light at a wavelength of λ = 520 nm, two lenses with focal lengths f = 150 mm, an iris aperture and a 2048 × 2048 CCD camera as a detector. The wing of a fly serves as a quasi-transparent object being placed at the distance r from the first lens where f < r < 2f. The iris diaphragm is positioned at the back focal plane of the first lens in order to limit the numerical aperture. Images were acquired at three different focal planes namely, z = Δz, z = 0 and z = −Δz where the defocus step is Δz = 1 mm. Figure 9(c) depicts the intensity variation along the optical axis which is computed from the under-focused (Fig. 9(b)) as well as the over-focused (Fig. 9(a)) image by the finite difference $d{I}_z^{\exp }=\frac{I(\mathrm{\Delta }z)-I(-\mathrm{\Delta }z)}{2\mathrm{\Delta }z}$.

 

Fig. 8 Schematic of optical setup.

Download Full Size | PDF

 

Fig. 9 Experimental images: a) Under-focused, b) Over-focused, c) Intensity variation along the optical axis.

Download Full Size | PDF

Figure 10(b) shows the reconstructed phase obtained by means of the GFA while Fig. 10(a) shows a Tikhonov regularization FFT-based solution which included the same padding used for the GFA but did not constrain the gradient of the phase. As clearly shown, the Tikhonov regularization reconstruction suffers from low-frequency artifacts due to the periodicity of the field of view as well as missing low frequency information. The red line across the retrieved phase highlights the fact that the FFT-based reconstruction experiences a steep slope in the empty area where the phase should be flat (Fig. 10(d)) while the gradient-constrained solution is reasonably constant in the absence of any object and also in places where we expect a constant wing thickness, as shown in Fig. 10(e). The TVAL3 solution shown in Fig. 10(c) fails to recover any low spatial frequency information within the wing at all.

 

Fig. 10 Reconstructed phase by a) Tikhonov regularization FFT-based, b) GFA method, c) TV minimization approach. (d),(e) and (f) Line profiles extracted along the red lines in each of the phase maps graphically depicted above.

Download Full Size | PDF

5. Conclusions

In this paper, we compared two algorithms for removing low spatial frequency artifacts in solutions of the TIE by constraining them to be sparse in either the ℓ₀ norm or the ℓ₁ norm of their gradient, namely, gradient flipping and total variation minimization as implemented by TVAL3. While the latter imposes convex constraints in the spirit of Compressed Sensing, the former solves a non-convexly constrained optimization problem. Both algorithms were able to recover piece-wise constant phase maps with and without periodic boundary conditions. In the case of noisy measurements, the GFA provided a solution which agreed slightly better with the input data at a high noise level of 30%, while performing slightly worse at a lower noise level of only 10%. For partially piece-wise constant and piece-wise linear test objects, and for Fresnel-propagation based measurements, the GFA solution agreed much better with the original phase used to simulate the test data. In a test on experimental data the GFA provided a physically very reasonable solution, while the TV-constrained solution did not seem physically reasonable.