1. INTRODUCTION

The imaging of samples at macro, micro, and nano scales is important in many fields of research, from physics (fluid mechanics, materials) to biology (cells, bacteria, viruses). Imaging remains challenging, and advanced techniques are still emerging. These new techniques often exploit physical principles involved in light–matter interactions from a variety of radiation sources (x-rays, electrons, visible or near-visible light, lasers). In this context, imaging techniques aim to record the perturbation of the incident electromagnetic wave by the sample of interest. The perturbed wave depends on both the absorption and phase-shift properties of the sample. Thus, being able to precisely measure the shape of this complex wave—the amplitude and phase of the electromagnetic field in space and time—is crucial to fully quantitatively characterize the sample.

Holography, invented by Gabor [1], is a technique to reconstruct the complex wave (amplitude and phase) due to the diffraction of light when a scattering sample is illuminated by a coherent source. This principle gave rise to a wide variety of imaging techniques (off-axis holography, phase shifting, in-line holography, diffractive tomography, x-ray diffractive microscopy, etc.) [2]. Since only intensity measurements of light are available, the problem of retrieving the phase on the sensor has fueled many studies. In contrast to interferometric setups such as off-axis or phase-shifting holography, iterative phase retrieval techniques are algorithms that estimate the phase of light in the plane of in-line holograms.

Digital hologram reconstruction is typically performed by backpropagating the hologram from the sensor plane to the object plane. In the absence of phase information with in-line holograms, the reconstruction suffers from artifacts called the twin image. To suppress the twin image, phase retrieval techniques were introduced in the 1970s and 1980s, and since then, improved algorithms have regularly been described in the literature. Most phase retrieval techniques are still derived from the methods of alternating projections initially proposed by Gerchberg and Saxton [3] and popularized and extended by Fienup [4,5]. This class of methods is still widely used today [6–13], with improvements to enforce a priori knowledge (support of the objects, admissible values domain, sparsity constraints) [14–16].

Inverse problems approaches take a different point of view from phase retrieval techniques: rather than recovering the phase on the sensor plane (which does not completely solve the sample reconstruction problem), they focus on the reconstruction of the complex-valued transmittance in the object plane. Because of measurement noise and hologram truncation at the borders of the sensor, backpropagating the complex amplitude in the sensor plane does not lead to perfect sample reconstruction, i.e., restoring the sensor phase and reconstructing the complex wave in the sample plane are not strictly equivalent problems.

Figure 1 gives an overview of the different ingredients of the inverse problems methodology in the context of holographic imaging. The “direct model” relates a description of the “parameters of interest” (i.e., the sample) and the “measurements” (i.e., the hologram), taking into account the “instrumental parameters” (i.e., the imaging setup). The sample may be described as a complex-valued transmittance plane [17], a volume of transmittance planes [18], or as a collection of objects [19]. Depending on the application, the sample is absorbing, and/or translucent (it induces a phase shift due to a refractive index difference with respect to the surrounding medium). It is described by the spatial distribution of the absorption, phase shift, or complex refractive index within a plane or a volume. When individual objects are considered, their 3D location, size, and complex refractive index are the parameters that characterize the objects. Any appropriate optical model can then be used to describe how an incident illumination wave gets diffracted by the sample and propagates to the sensor: diffraction by a plane (Rayleigh–Sommerfeld, Fresnel approximation, angular spectrum) [20], diffraction by a volume (Born approximation, Rytov approximation) [21], diffraction by a sphere (Mie scattering) [22], or free-space propagation. Multiple holograms can be recorded to improve information diversity by including lateral or axial shifts of the sensor [16,23], multiple illumination wavelengths [24], or multiple illumination angles [25,26]. When a suitable model of the observations has been defined, a “reconstruction” (or “estimation”) method can be selected. If the sample is described by a spatial distribution within a plane or a volume, the problem amounts to an “image or volume reconstruction.” If a collection of objects is to be identified, the geometrical and optical parameters of each object have to be estimated (“parameter estimation”). “Image or volume reconstruction” problems are generally solved by regularized inversion [27–37], i.e., by solving a minimization problem whose objective function is defined by a data-fitting term that penalizes the deviation of the model from the data. Regularization terms are added to favor reconstructions with desirable properties (e.g., smooth images, sharp edges, zero-valued background). Additional (hard) constraints are also enforced to guarantee that the reconstruction fulfills basic physical properties such as the positivity of the absorption or the sign of the phase shift induced by the objects (e.g., thin objects with an optical index larger than that of the surrounding medium induce a positive phase shift). Depending on the type of cost function to minimize [quadratic, nonlinear but differentiable (smooth), non-differentiable (non-smooth), non-convex], different strategies can be used. “Parameter estimation” of a collection of objects generally follows a different path because the reduced number of unknowns makes the regularization unnecessary. Rather than considering a fixed discretization of space, the 3D location of each object is generally sought in a continuous domain (i.e., subpixel location). Objects are often detected and characterized one after another, and the measurements are described by the superimposition of the diffraction patterns due to each object [19,38]. Nonlinear least-squares minimization techniques are applied to estimate the geometrical and optical parameters of each object. In the case where the data contain some signatures that cannot be fitted by the model (outliers), e.g., out-of-field objects, the estimations can be improved by using a robust signal processing approach such as the reweighted least squares [39].

 

Fig. 1. Hologram reconstruction based on inverse problems approaches: the direct model connects the sample of interest to the measurements (the holograms); inverting this model leads to a reconstruction and/or estimation of the optical parameters of the sample. The case chosen as an illustration in this paper is highlighted in yellow.

Download Full Size | PDF

The hologram formation model depends on “instumental parameters” such as the wavelength and coherence of the source, the sample-to-sensor distance, or the sensor response (fill-factor and spectral sensitivity). The model of the data not only underlies the reconstruction algorithm, but also provides a way to calibrate the experimental parameters.

The goal of this paper is to illustrate the inverse problems methodology and to draw connections with the popular phase retrieval techniques based on Fienup’s approach. We selected a problem that can be addressed using either approach: the reconstruction of an in-line hologram of a non-absorbing plane sample that induces spatially varying phase shifts. We discuss the advantages provided by the inverse problems framework and interpret Fienup’s as an instance of a proximal gradient algorithm applied to a specific cost function.

The paper is organized as follows: next, we present the phase retrieval problem in the context of in-line digital holography. We then present the standard alternating projections strategy proposed by Gerchberg and Saxton, Fienup, and subsequent variants. In Section 4, we present the inverse problems methodology applied to the reconstruction of in-line holograms of phase-only samples. We validate and compare the algorithms in Section 5, first in numerical simulations, then on experimental in-line holograms of microscopic beads. Finally, we discuss the possible extensions of inverse problems approaches for phase retrieval.

2. THE INVERSE PROBLEM OF IMAGE RECONSTRUCTION IN IN-LINE HOLOGRAPHY

The object of interest, i.e., the sample, can be modeled as a two-dimensional (2D) complex transmittance plane $ \underline t ({\boldsymbol r}) $ [$ {\boldsymbol r} = (x,y) $ is the vector of 2D spatial coordinates, and all complex-valued variables are underlined in this paper]. A transmittance $ \underline t ({\boldsymbol r}) = 1 $, $ \forall \,\,{\boldsymbol r} $ means a fully transparent plane that leaves any illumination wave unperturbed after passing through it. Consequently, the presence of scattering objects in this plane will be revealed by a certain deviation $ \underline o ({\boldsymbol r}) $ from a unit transmittance:

Digital in-line holography (DIH) (cf. Fig. 2) consists of illuminating a sample plane with a coherent (generally plane) wave $ {\underline a _0}({\boldsymbol r}) $ of wavelength $ \lambda $, and recording on a digital sensor the intensity $ I({\boldsymbol r}) $ of the total diffracted wave $ {\underline a _z}({\boldsymbol r}) $ at a distance $ z $ from the sample plane (the sensor plane being parallel to the sample plane). From the Huygens–Fresnel principle, the physical data acquisition process is written as

 

Fig. 2. In-line holography principle: a plane wave illuminates the sample plane with normal incidence. The presence of inhomogeneity in the spatial distribution of the complex-valued refractive index, in the sample plane, distorts the illumination wave. A diffraction pattern is recorded on the sensor plane.

Download Full Size | PDF

The backpropagation kernel $ {\underline h _{ - z}} $ can be used to recover the original complex wave from a diffracted (i.e., propagated) wave.

Considering Eq. (1), and under the assumption that the incident wave is plane ($ {\underline a _0}({\boldsymbol r}) = {\underline a _0} $), Eq. (2) becomes

The intensity $ I({\boldsymbol r}) $ is recorded by a digital sensor. Therefore, effective intensity measurements correspond to sampled (and noisy) intensity values that can be collected in the form of a column vector $ {\boldsymbol d} = ( {{d_q}} )_{q = 1 \ldots M}^{T} $ with $ M $ data measurements. This also leads to considering the targeted physical quantity $ \underline o $ as a vector $ \underline {\boldsymbol o} = ( {{x_q}} )_{q = 1 \ldots N}^{^{T}} $ of $ N $ unknowns. In imaging systems, $ N $ is often equal to $ M $, as we aim to recover an image of the same size as the hologram, i.e., the field of view. In the context of inverse approaches, the unknown field of view of $ \underline {\boldsymbol o} $ can be extended, leading to $ N \gt M $ (cf. Section 6). We denote $ {\mathbb D} $ and $ {\mathbb O} $, respectively, as the feasible domains of $ {\boldsymbol d} $ and $ \underline {\boldsymbol o} $, i.e., the domains of valid values for each quantity. From Eq. (4), $ {\mathbb D} \subset {{\mathbb R}^M} $ and $ {\mathbb O} \subset {{\mathbb C}^N} $. As discussed in Section 1, in order to ensure physically relevant properties of the object of interest, the feasible domain $ {\mathbb O} $ can be restricted, for example, to impose a positive absorption or the sign of the phase shift. Such properties are applied in the regularized inverse approach developed in our case study in Section 4.B.

The algebraic process that “transforms” $ \underline {\boldsymbol o} $ into $ {\boldsymbol d} $ is modeled by the mapping $ {\boldsymbol m} $, which goes from $ {\mathbb D} $ to $ {\mathbb O} $, and is a discretized version of the physical—analytic—formulation $ m $ in Eq. (4). Finally, we get the following algebraic formulation of the image formation model:

Retrieving the unknown image of the sample $ \underline{\boldsymbol o} $, denoted $ {\underline{\boldsymbol o} ^*} $, from the observed hologram $ {\boldsymbol d} $, based on the hologram formation model in Eq. (5), constitutes the inverse problem of holographic reconstruction.

3. FIENUP’S ALTERNATING PROJECTIONS STRATEGIES

The standard alternating projections strategy, still widely used to date, was first proposed by Fienup [4] as the error-reduction (ER) algorithm. It is an upgraded variant of the initial method proposed by Gerchberg and Saxton [3] from two-plane intensity measurements.

In Ref. [5], Fienup has proposed several variants of the ER algorithm, such as the basic input–output (BIO) or the hybrid input–output (HIO), involving a parameter $ \beta $ in the “projection on data” step that relaxes the strict projection on the feasible domain. These variants are known to be more efficient than the classical ER strategy (HIO is considered as the most efficient). In this paper, we focus on the classical approach ER, which we find the most intuitive way to understand the alternating projections strategy. We show in Section 4.C.3 that this algorithm can be reformulated to fit the inverse problems framework. Thus, our goal is not an exhaustive comparison of inverse approaches with all variants of the alternating projections strategy, but rather to give an “inverse problems” interpretation of Fienup’s strategy and show that inverse problems offer a framework that is accessible, yet more flexible and powerful than alternating projections.

Algorithm 1 summarizes the Fienup ER strategy. Details of implementation are given in Appendix A, in Algorithm 3.

Fienup’s ER algorithm requires only one intensity measurement image $ {\boldsymbol d} $ at the sensor plane at distance $ z $. This measurement image has to be normalized so that the background equals 1, to ensure that the retrieved transmittance plane $ \underline{\boldsymbol o} $ is normalized the same way ($ \underline{\boldsymbol o} = \textbf{1} $ stands for a fully transparent plane; see Section 2), and that appropriate physical constraints can be applied. The normalized data image, denoted $ {\boldsymbol{ \bar d}} $, can be obtained, for example, by dividing $ {\boldsymbol d} $ by its mean or its median (almost valid in the case of a low-density distribution of objects).

Algorithm 1. Fienup’s ER algorithm [4]

View Table | View all tables in this article

The principle is to iteratively alternate the following steps starting from an estimate of the object $ {\underline{\boldsymbol o} ^{(i)}} $ at the sample plane to get a new estimate $ {\underline{\boldsymbol o} ^{(i + 1)}} $:

This approach is based on the approximation $ \underline{\textbf{H}} _z^\dagger {\underline{\textbf{H}}_z} \approx \textbf{I} $ that makes it possible to go back and forth between the sample plane and the sensor plane. Since $ \underline{\textbf{H}} _z^\dagger $ is not the exact inverse of $ {\underline{\textbf{H}} _z} $, there are strong limitations related to sampling and truncation effects (e.g., objects on the border of the field of view are inevitably distorted). If the algorithm is stopped after several iterations and the sensor plane constraint applied, the data are augmented by the phase information (i.e., the lost phase on the sensor plane is retrieved by the algorithm). Due to this restored phase, the twin-image phenomenon typical of in-line holography is strongly reduced.

The great popularity of this class of algorithms is due to the fact that the strategy is very intuitive as well as being easy to implement. Likewise, these algorithms have been proven to converge to a local minimum [40–42]. Such approaches have been exploited mainly in the fields of x-ray coherent diffraction imaging [7,43–46] and digital holographic microscopy [29,47,48].

As a result, research in this field is still active and new approaches are often proposed. In the past 3 decades, it has benefited from many practical improvements [6,10,11,13,49] and theoretical studies [9,41,50]. Recent approaches include new constraint enforcement strategies using sparsity constraints in the object or data spaces [51,52] or in the wavelets domain [16,53].

4. INVERSE PROBLEMS METHODOLOGY FOR THE RECONSTRUCTION OF IN-LINE HOLOGRAMS: A TUTORIAL

A. Linearized Image Formation Model

The first step of the inverse problems methodology is to derive a suitable image formation model, by analyzing the object and data characteristics. Here, we describe a linearized model based on two hypotheses: (i) the sample is a distribution of purely and weakly dephasing objects [$ \rho ({\boldsymbol r}) = 1 $; see Eq. (1)], such that

B. Building a Regularized Inverse Approach

First, we recall the definition of the $ {l_2} $ and $ {l_1} $ norms applied to a vector $ {\boldsymbol x} $: $ ||{\boldsymbol x}{||_2} = \sqrt {\sum\nolimits_q x_q^2} $ and $ ||{\boldsymbol x}{||_1} = \sum\nolimits_q |{x_q}| $. The notation $ ||{\boldsymbol x}{||_{\textbf{W}}} $, where $ \textbf{W} $ is a positive definite matrix, corresponds to $ ||{\boldsymbol x}||_{\textbf{W}}^2 = {{\boldsymbol x}^{T}}\textbf{W}{\boldsymbol x} $.

1. Data Fidelity

An inverse approach consists of retrieving an optimal solution $ {{\boldsymbol o}^*} $ to Eq. (8), knowing an approximate numerical model $ {\boldsymbol{\tilde m}}({\boldsymbol o}) $ of the data formation process.

To do so, we define a calculable metric $ {{\cal J}_{\textsf{fid}}}({\boldsymbol o},{\boldsymbol d}) $ that evaluates the error between the data $ {\boldsymbol d} $ and the model $ {\boldsymbol{ \tilde m}}({\boldsymbol o}) $ obtained from a guessed $ {\boldsymbol o} $. This error term $ {{\cal J}_{\textsf{fid}}}({\boldsymbol o},{\boldsymbol d}) $, also called the data-fidelity term, acts as a penalization: the larger the error, the farther the guess $ {\boldsymbol o} $ is from the solution and the more it has to be corrected to reach a better estimate.

The most popular data-fidelity term is the weighted least-squares criterion, which derives from the assumption of Gaussian errors. The criterion corresponds to the $ \textbf{W} $-weighted squared $ {l_2} $ norm of the differences between $ {\boldsymbol d} $ and $ {\boldsymbol{ \tilde m}}({\boldsymbol o}) $, which is written as

To solve the inverse problem, the objective is to find the solution $ {{\boldsymbol o}^*} $ that minimizes the criterion Eq. (9), thus the error between the data $ {\boldsymbol d} $ and the model $ {\boldsymbol{\tilde m}}({{\boldsymbol o}^*}) $. Such a minimization problem is written as

2. Constraints and Regularizations

The minimization problem Eq. (10) is not sufficient to obtain a satisfactory solution because measurements are noisy, and there are too many unknowns for a reliable estimation based only on data fitting. As a result, unsatisfactory reconstructions can be found that perfectly match the data, i.e., solutions of the problem where the noise is also fitted by the model.

These issues are the characteristics of what is called an ill-posed problem. To overcome these limitations, the resolution of the problem cannot be achieved by considering only the data-fidelity term Eq. (10). One has to find a solution $ {{\boldsymbol o}^*} $ that best fits a particular prerequisite on the targeted information: one has to enforce some prior knowledge about the unknown $ {\boldsymbol o} $.

Such prior knowledge can be injected in the minimization problem Eq. (10), taking the form of hard and/or soft constraints. The new minimization problem is written as

$ {\boldsymbol o} \in {\mathbb O} $ imposes the feasible domain of values $ {\mathbb O} $ for the estimate $ {\boldsymbol o} $. In many inverse problems involving real quantities, a relevant constraint is to enforce positivity. This is, for example, the case in positron emission tomography (PET), where the quantity of interest is a photon count [57].

$ {{\cal J}_{\textsf{reg}}}({\boldsymbol o},{\boldsymbol \theta }) $ are soft constraints imposed as a combination of regularization terms, parameterized by the set of hyperparameters $ {\boldsymbol \theta } $ that balance each term with respect to the data-fidelity term. In general, standard regularizations aim to smooth the estimate $ {\boldsymbol o} $ to filter high frequencies due to noise. In the following, we present three popular regularizations in image processing. Figure 3 illustrates their behavior on reconstructions, using a regularized inverse problems approach, of the simulated hologram that will be used in Section 5 for further demonstrations. More details on the method, the simulation, and reconstruction parameters are also provided in Section 5.

 

Fig. 3. Illustration of the behavior of several regularization terms on reconstructions (in a region of interest) of a simulated hologram (a) using a regularized inverse problems approach. Details on the method, the simulation, and reconstruction parameters are given in Section 5 and Tables 1 and 2. (a) Data (in-line hologram). The yellow frame indicates the region of interest that is extracted from the field of view for visualization. (b) Ground truth phase. (c)–(f) Reconstructed images using the following regularizations: (c) squared $ {l_2} $ norm of the gradient $ {{\cal J}_{{{l}_{2}}\nabla }} $ [Eq. (15)], (d) edge-preserving smoothness $ {{\cal J}_{{{\text{TV}}_ \epsilon }}} $ [Eq. (16)], (e) $ {l_1} $ sparsity constraint $ {{\cal J}_{{l_1}}} $ [Eq. (17)], and (f) composition of the $ {l_1} $ sparsity constraint and edge-preserving smoothness $ {{\cal J}_{\textsf{reg}}} = {{\cal J}_{{l_1}}} + {{\cal J}_{{{\text{TV}}_ \epsilon }}} $ [Eq. (18)]. Red curves show a line profile passing through two particles of the image.

Download Full Size | PDF

A possible strategy is to minimize the squared $ {l_2} $ norm of the gradient image $ \nabla {\boldsymbol o} $:

Direct sparsity of the image can also be enforced. This strategy is very often used to favor a low-density distribution of objects in the image. This is implemented by minimizing the separable $ {l_1} $ norm of the image:

These regularizations can be mixed to benefit from each of their advantages [cf. Fig. 3(f)], leading $ {{\cal J}_{\textsf{reg}}}({\boldsymbol o},{\boldsymbol \theta }) $ in Eq. (14) to be a combination of some of these terms. This increased flexibility comes at the cost of a more difficult tuning of the reconstruction algorithms because of the increase in the number of hyperparameters. For example, in Fig. 3(f) and reconstruction experiments in Section 5, we use a combination of a separable sparsity prior $ {{\cal J}_{{l_1}}} $ [Eq. (17)] and an edge-preserving regularization $ {{\cal J}_{{{\text{TV}}_ \epsilon }}} $ [Eq. (16)], leading to the following expression for the global regularization term $ {{\cal J}_{\textsf{reg}}}({\boldsymbol o},{\boldsymbol \theta }) $:

C. Case Study: Solving a Standard Sparsity-Based Problem

1. Problem

Following the methodology presented in the previous section, we now focus on the resolution of the following criterion that corresponds to a popular sparsity-based approach:

In the following, we introduce a simple yet efficient method to solve this inverse problem, named Iterative Shrinkage-Thresholding Algorithm (ISTA). This method belongs to the class of proximal gradient methods [60], whose general principle is depicted in Fig. 4. The key idea is to replace the original minimization problem by a sequence of simpler minimization problems involving a surrogate function that majorizes the original function. The cost function is decomposed into the sum of a smooth (i.e., differentiable) and a non-smooth (i.e., non-differentiable) term. The smooth term is replaced by a separable quadratic approximation. Minimizing this approximation corresponds to computing the proximal operator associated with the non-smooth term [60]. In the case of regularization by an $ {l_1} $ norm, under a positivity constraint, we recall in Fig. 5 that the proximal operator corresponds to a soft-thresholding operation with a clipping of negative values to zero. In the following, we also establish that there is a strong connection between Fienup’s alternating projections method and this strategy applied to a particular formulation of the minimization problem.

 

Fig. 4. Principle that underlies proximal gradient methods is the iterative minimization of a local approximation of the original cost function. If the parameter $ t $ is chosen small enough, the approximation is majorant, and minimizing the approximation improves the current solution until convergence. When the non-smooth component $ {\cal H} $ is separable (as is the case of the $ {\ell _1} $ norm), minimizing the local approximation is easily done [see Fig. (5)].

Download Full Size | PDF

 

Fig. 5. Computation of the proximal operator: an illustration with the case of the $ {l_1} $ norm under positivity constraints. Values that are smaller than $ \mu t $ are mapped to zero by the proximal operator. Larger values of the regularization weight $ \mu $ therefore lead to more values being set to zero, i.e., a sparser solution.

Download Full Size | PDF

 

Fig. 6. Phase image reconstructed from simulated data with different methods and regularization terms. (a) Data (in-line hologram). The yellow frames indicate the regions of interest that are extracted from the field of view for visualization. (b) Ground truth phase. (c) Evolution of the normalized cost criterion (minimum and maximum values ranged between zero and one) for each reconstruction. (d)–(i) Reconstructed images using (d) Fienup method, (e) Fienup method + soft-thresholding, (f) Fienup method + soft-thresholding, stopped at 10 iterations, (g) RI method using FISTA with soft-thresholding, (h) RI method using FISTA with soft-thresholding + edge-preserving, (i) RI method using FISTA with soft-thresholding + edge-preserving, stopped at 10 iterations. For each reconstruction, a positivity constraint is imposed on the solution. The values of the hyperparameters are indicated in Table 2. Black dots delimit the contour of the ground truth image. Red curves show line profiles passing through particles of the image.

Download Full Size | PDF

 

Fig. 7. Phase image reconstructed from experimental data with different methods and regularization terms. (a) Data (in-line hologram). The yellow frame indicates the region of interest that is extracted from the field of view for visualization. (b) Evolution of the normalized cost criterion (minimum and maximum values ranged between zero and one) for each reconstruction. (c)–(h) Reconstructed images using (c) Fienup method, (d) Fienup method + soft-thresholding, (e) Fienup method + soft-thresholding, stopped at 10 iterations, (f) RI method using FISTA with soft-thresholding, (g) RI method using FISTA with soft-thresholding + edge-preserving, (h) RI method using FISTA with soft-thresholding + edge-preserving, stopped at 10 iterations. For each reconstruction, a positivity constraint is imposed on the solution. The values of the hyperparameters are indicated in Table 2. Red curves show a line profile passing through two particles of the image.

Download Full Size | PDF

2. Optimization Strategy

The application of the proximal gradient framework to the minimization of Eq. (19) leads to the following iteration:

To ensure the convergence of the algorithm, the step length $ t $ has to be chosen adequately. For our problem, a suitable value depends on the maximum eigenvalue of $ \textbf{G}_z^{T}{\textbf{G}_z} $ (i.e., the Lipschitz constant of $ {{\cal J}_{\textsf{fid}}} $; see [60,61]). As $ {\textbf{G}_z} $ is a convolution operator, this is equivalent to the maximum squared modulus of the Fourier transform of the kernel $ {g_z} $. It is also possible to use back-tracking methods to reduce the step length $ t $ if it is found that the surrogate function does not majorize the original cost function for the current value of $ t $ [61].

In Ref. [61], the authors propose a strategy to accelerate the convergence rate of ISTA, leading to the popular algorithm named Fast Iterative Shrinkage-Thresholding Algorithm (FISTA). We detail the general steps of FISTA in Algorithm 2. Implementation details are given in Appendix A, Algorithm 4. In Algorithm 2, if the scalar factor $ s $ is kept to the value 1, we fall back to the simpler ISTA algorithm.

Algorithm 2. FISTA algorithm [61] for solving problem Eq. (19)

View Table | View all tables in this article

3. Alternating Projections Point of View

We analyze the alternating projections strategy under the ISTA formulation presented in Fig. 4. From steps 1 to 4 in Algorithm 1, we derive the equivalent iteration step of Fienup’s algorithm as follows:

The term in the square brackets on which the projector $ {{\cal P}_{\mathbb O}} $ is applied corresponds to a steepest gradient descent step that decreases the following data-fidelity criterion:

The projector $ {{\cal P}_{\mathbb O}} $ is the orthogonal projection onto the convex set $ {\mathbb O} $, and corresponds in our case to forcing positive or negative values of the real and/or imaginary parts of $ \underline {\boldsymbol o} $, and/or to forcing to zero all pixels that are outside the support of the object. $ {{\cal P}_{\mathbb O}} $ is in fact the proximal operator associated with the indicator function $ {\iota _{\mathbb O}} $ of the domain $ {\mathbb O} $ [60]:

Then, applying Fienup’s alternating projections strategy given in Eq. (23) is analogous to performing ISTA iterations with a fixed step length $ t = \frac{1}{2} $, derived from the following problem:

The above analysis is important since it demonstrates that the alternating projections method underlies a particular inverse problems approach, and thus falls within this rigorous framework. The same as with the method developed in Section 4.C.2, the accelerated FISTA algorithm can then be applied to this problem, and it is even possible to replace the standard projection $ {{\cal P}_{\mathbb O}} $ by, for example, a soft-thresholding proximal operator $ {{\cal T}_{\mu t}} $ to enforce sparsity of $ \underline {\boldsymbol o} $ (real or/and imaginary part), or other regularization terms (provided that an efficient computation of the associated proximal operator is possible).

This interpretation of the alternating projections strategy in the framework of inverse problems and convex optimization methods has already been noted in previous works. In Ref. [5], Fienup already analyzed his variants of the ER algorithm as particular steepest descent strategies. Levi and Stark [50] have established that this method can be analyzed as a nonconvex instance of the Projection Onto Convex Sets (POCS) algorithm. This was rigorously demonstrated by Bauschke et al. [41], since Fienup’s BIO and HIO are shown to be instances of Dykstra and Douglas–Rachford algorithms. In Ref. [12], alternating projections strategies are reformulated as proximity operators derived from the maximum likelihood point of view.

Conversely, this analysis makes it possible to interpret the ISTA strategy in Eq. (20) from the alternating projections point of view, since we can deduce that this algorithm performs a succession of a propagation step $ {\textbf{G}_z} $ of the image $ {\boldsymbol o} $ to get a model of the data $ {\boldsymbol{ \tilde m}}({\boldsymbol o}) $, followed by a pseudo-backpropagation step $ \textbf{G}_z^{\text{T}} $ of an error term [difference between $ {\boldsymbol{ \tilde m}}({\boldsymbol o}) $ and $ {\boldsymbol d} $] and finally a projection step $ {{\cal T}_{\mu t}} $ on constraints.

Table 1. Experimental and Simulation Parameters of In-Line Holograms Data

View Table | View all tables in this article

Table 2. Reconstruction Parameters for Simulations and Experiments in Figs. 3,6–8

View Table | View all tables in this article

 

Fig. 8. Reconstruction (detection) of an out-of-field particle using RI method and field-of-view extension. (a) Data in-line hologram highlighting (in yellow) a part of the out-of-field hologram. (b) Ground truth image. (c), (d) Reconstructions with two different sets of regularization weights values [cf. Table 2]. The red frames show the initial data field of view (pixels seen by the detector). The green frames are zooms on the out-of-field particle.

Download Full Size | PDF

With the above interpretation, we want to demonstrate that the inverse problems framework, in some specific but adapted cases, can keep the intuitive nature and ease of implementation property that generally constitutes the argument for preferring the use of Fienup’s method. This is an original contribution of this paper compared with the previous works.

5. RESULTS

In this section, we present the reconstruction results obtained with the proposed inverse problems approach formulated in Eq. (19), using the ISTA iteration given in Eq. (20) and its acceleration with the FISTA algorithm Algorithm 2. We compare this method with Fienup’s alternating projections approach defined by the iteration Eq. (23) using two versions: (i) the standard way with positivity constraint enforcement and (ii) an upgraded way with a soft-thresholding step ($ {l_1} $ sparsity constraint). All the reconstruction strategies also enforce a positivity constraint. In the following interpretations of the results, our regularized inversion approach method is referred to as “the RI method,” while Fienup’s alternating projections approach is referred to as “the Fienup method.”

 

Fig. 9. Derivation of a reconstruction algorithm from a cost function.

Download Full Size | PDF

A. Data

We reconstruct two data sets: an experimental in-line hologram of polystyrene beads in immersion oil [cf. Fig. 7(a)], and a simulated in-line hologram of beads in similar conditions [cf. Fig. 6(a)]. In this simulation, the transmittance $ {\boldsymbol o} $ is modeled as 2D truncated Gaussian footprints (the Gaussians are truncated so that the footprint has a diameter of $ 6\sigma $, and zero outside.). The peak value of this footprint is set to a refractive index difference $ ({n_0} - {n_{\text{bead}}}) = 0.07 $. The in-line hologram is simulated under Fresnel approximation. Table 1 summarizes the experiment and simulation parameters.

Looking at the data in Fig. 6(a), we can see that a particle outside the field of view is also present. Its diffraction fringes are truncated at the top of the hologram. We intentionally included this out-of-field object to illustrate that the inverse approach framework allows the reconstruction of holograms in a wider field of view than that seen by the sensor. We discuss this particular point in Section 6.

In these conditions, the maximum phase shift induced by the particles is equal to $ \Delta \varphi = 2\pi {d_{\text{bead}}}({n_0} - {n_{\text{bead}}})/\lambda \approx 0.83\,\,\text{rad} $. Then, the hypotheses (weakly dephasing objects) made to derive the model proposed in Section 4.A are not valid. Still, in addition to the fact that the problem is simplified (only a real image has to be estimated), we show that this approximate model leads to good reconstructions.

B. Reconstructions

For each experiment, all the reconstruction parameters in Figs. 6 and 7 are summarized in Table 2.

Figure 6 shows the results of the reconstruction of the simulated hologram data, using the two methods presented above, with different parameters. For each reconstructed image in Figs. 6(d)–6(i), the values of hyperparameters are given in Table 2. First, we can observe that both reconstruction methods using a soft-thresholding step give a satisfying estimation of the objects’ support. We see that the Fienup method enforcing this sparsity constraint is clearly better than the standard approach with a positivity constraint alone. However, the simulated objects have a spatially extended support, i.e., they cannot be considered as point objects. Thus, the sparsity constraint tends to create tiny holes within the objects. The sparsity constraint alone is not well adapted to this particular hologram reconstruction problem. The RI method can be improved by adding an edge-preserving regularization in the form of the edge-preserving term presented in Eq. (16). The new minimization problem is written as

The reconstruction via this combination of regularizations is shown in Fig. 6(h), where a clear benefit of considering both a sparsity constraint and an edge-preserving smoothing constraint is observed, as already shown in Fig. 3(f). Using the inverse problems interpretation of the Fienup method, we could also adapt Fienup’s algorithm to include both regularizations.

In Fig. 6(c), we show the evolution over the iterations of the global (normalized) criterion of the RI and Fienup methods, respectively, the criteria Eqs. (27) and (26). We illustrate the convergence of the RI method with several optimization strategies: FISTA, ISTA, and VMLMB [64,65], a quasi-Newton gradient variable metric method with limited memory requirements and possibly bound constraints enforcement on the unknowns. Whatever the method, an empirical convergence of the criteria values is reached in less than 10 iterations. This is confirmed in the reconstructions Figs. 6(f) and 6(i), which show the same reconstructions as Figs. 6(e) and 6(h) after only 10 iterations. We see that both reconstructions (RI and Fienup) are almost the same as the solution obtained after 1000 iterations.

Figure 7 shows reconstructions of the experimental hologram data, using the two methods presented above with different parameters. The field of view has been cropped to $ 1080 \times 1080 $ images. For each reconstructed image in Figs. 7(c)–7(h), the values of the hyperparameters are given in Table 2. In Fig. 7(b), the evolution over the iterations of the global normalized criterion of the methods is shown.

The observations and analysis that can be performed on these reconstructions match those made on the simulated data: the objects of interest are clearly detected. Moreover, the estimations of their shape, diameter [see Fig. 7(h)], and phase shift are coherent with the expected values. Indeed, from the experimental parameters in Table 1, the maximum phase shift should be around 0.83 rad. Again we notice that the convergence of the criteria is almost reached in a few tens of iterations.

6. DISCUSSION

As pointed out in Section 5.A, the linear model derived in this paper is only an approximation. Nevertheless, we observe in both the simulated and experimental results that the shape of the beads is correctly retrieved and that the phase-shift values estimated remain close to the expected values (the difference is less than one order of magnitude). The values obtained on the experimental holograms even match the expected phase shift at the center of the beads ($ \Delta \varphi \approx 0.83\,\,\text{rad} $).

This illustrates that the regularizations and constraints introduced to solve the inverse problem balance the modeling errors and lead to a satisfying solution: the reconstruction is robust to measurement and modeling errors.

In Section 5.A, we mentioned that an out-of-field particle was considered in the simulation, leading to cropped diffraction fringes in the upper part of the data [cf. Fig. 8(a)]. The inverse problems framework makes it possible to account for this data truncation (a problem also frequently faced in x-ray computerized tomography [66]) and reconstruct a wider field of view. This can be done by two equivalent ways: (i) an extension of the size of both the reconstructed object plane and the measured hologram; (ii) a field extension of only the reconstructed object plane. In the first case, unmeasured pixels (i.e., pixels outside the field of view of the sensor) are given a weight $ {w_q} = 0 $ in Eq. (9). In the second case, the model is rewritten to include a truncation operator (i.e., a rectangular matrix obtained by chopping off rows corresponding to unmeasured pixels):

Algorithm 3. Detailed implementation of Algorithm 1

View Table | View all tables in this article

Algorithm 4. Detailed implementation of Algorithm 2

View Table | View all tables in this article

Figure 8 shows reconstructions of the simulated hologram with the RI method, where the object field of view is doubled compared to the data field of view. These reconstructions illustrate that the out-of-field particle can be detected. Since a major part of the diffraction pattern of this bead is missing, the reconstruction is not as good as for other beads. If a proper regularization is used, the bead can be detected (a “too large” sparsity regularization would tend to suppress it). Here, the signal-to-noise ratio and the approximate direct model lead to an imperfect reconstruction of the morphology of the bead (asymmetrical shape). Other “reconstruction” approaches based on a parametric model of the spherical beads coupled with an adequate inversion algorithm (e.g., continuous matching pursuit algorithm) lead to accurate estimations even out of the field of view [19].

These discussions give an overview of the refinements made possible by a properly built inverse approach to enhance the quality of reconstructions and extract the most relevant information from the data.

7. CONCLUSION

In this work, we have presented the inverse problems methodology applied to the reconstruction of the phase information from in-line intensity holograms. The main goal was to provide a tutorial for this reconstruction strategy, in comparison with the standard alternating projections method proposed by Fienup. To this end, we have presented the overall methodology for building an inverse approach dedicated to the targeted application, from the modeling of data formation to the choice of suitable constraints and regularizations. Deriving a reconstruction algorithm from a cost function is straightforward, for classical regularization terms, by following a decomposition into smooth and non-smooth components, as summarized in Fig. 9. We have shown that Fienup’s method is analogous to a particular formulation of this phase retrieval problem when solved by proximal gradient descent iterations, i.e., the ISTA algorithm. We have shown reconstructions from our proposed RI method compared with standard and upgraded versions of Fienup’s method, and for which we have provided algorithmic details of implementation. The analysis of the results has shown that the reconstruction quality can be quite similar when the Fienup method is refined using the inverse problems framework. We hope that we have convinced the reader that developing a reconstruction strategy following the inverse problems methodology can bridge a gap in reachable reconstruction quality, while allowing much more flexibility to extract relevant information from the data, at no cost in terms of implementation efforts and computational burden. We think that this unifying point of view on various approaches to hologram reconstruction will help the cross-fertilization of algorithmic ideas.