## Abstract

Reconstruction of phase objects is a central problem in digital holography, whose various applications include microscopy, biomedical imaging, and fluid mechanics. Starting from a single in-line hologram, there is no direct way to recover the phase of the diffracted wave in the hologram plane. The reconstruction of absorbing and phase objects therefore requires the inversion of the non-linear hologram formation model. We propose a regularized reconstruction method that includes several physically-grounded constraints such as bounds on transmittance values, maximum/minimum phase, spatial smoothness or the absence of any object in parts of the field of view. To solve the non-convex and non-smooth optimization problem induced by our modeling, a variable splitting strategy is applied and the closed-form solution of the sub-problem (the so-called proximal operator) is derived. The resulting algorithm is efficient and is shown to lead to quantitative phase estimation on reconstructions of accurate simulations of in-line holograms based on the Mie theory. As our approach is adaptable to several in-line digital holography configurations, we present and discuss the promising results of reconstructions from experimental in-line holograms obtained in two different applications: the tracking of an evaporating droplet (size ∼ 100μm) and the microscopic imaging of bacteria (size ∼ 1μm).

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

The principle of in-line holography was first proposed by Dennis Gabor in 1948 [1]. It uses a coherent - or semi coherent - optical illumination to image absorbing and/or phase objects. This is a method of choice for the imaging of transparent objects, which is impossible with conventional imaging.

While interferometric setups like off-axis holography [2] or phase shifting [3] can be used to measure the phase in the hologram plane and then recover the complex amplitude in the object plane by simple back-propagation of the complex wave field, with in-line holography a numerical inversion is necessary to recover the phase either in the hologram plane or in the object plane. Phase retrieval is a central problem in in-line holography [4].

State-of-the-art methods propose for example variants of the widely used Gerchberg-Saxton algorithm [5], such as the Error-Reduction Fienup’s algorithm [6] which includes a step for constraining the support of the object of interest in the object plane.

In the past decade, more sophisticated methods, using additional processing inspired by advances in the domain of compressive sensing, have been proposed, *e.g.* enforcing sparsity constraints in the object or data spaces [7,8] or in the wavelets domain [9–11]. Some methods use several holograms acquired at different distances to improve the amount of information [9,11]. Very recently, approaches based on deep learning have adapted image restoration techniques to phase retrieval [12,13].

Phase retrieval requires the injection of constraints on the holographic objects, mostly in the form of support constraints (objects of finite size, separated by empty regions). Correctly estimating the support of the objects would require a good pre-estimation of the phase, but the phase cannot be recovered without knowledge of the support: a chicken and egg problem that can only be addressed by jointly recovering the phase information and the objects support. Inverse problems approaches based on a fine model of hologram formation and including several constraints on the objects (based on application-specific knowledge) are well suited to perform this joint estimation. Several inverse problems methods have been proposed in digital holography [14–18] and applied to various phase retrieval problems : in-line holography [17,19,20], off-axis holography [14,16,18] and diffractive tomography [21]. A key point of such approaches is that they often need to solve a non-smooth and non-convex optimization problem. In this context, proximal algorithms can be applied to both convex [22,23] and non-convex problems [17,24].

The work presented here belongs to the latter. We propose a regularized inverse problems approach to reconstruct the real and imaginary parts of the object’s complex transmittance (i.e. phase and amplitude) starting from a single in-line hologram. In addition to using a non-linear direct model of image formation, our approach includes several physical constraints (*i.e. prior* information) on the transmittance plane, in the form of hard constraints (strictly enforced) and soft constraints (penalization terms to favor regularity or sparsity). Our optimization problem therefore includes non-linear constraints, smooth, non-smooth or non-convex terms. To tackle the difficulty involved in minimization, we use a “variable splitting” strategy and derive the closed-form global solution of the sub-problem that contains all the non-linear constraints and non-smooth terms, and part of the non-convexity. This proximal operator is an important contribution of our work since it forms the core of the reconstruction method and reduces the number of additional parameters involved in the variable-splitting process [25].

In section 2, we present the inverse problem for the reconstruction of absorbing and phase objects, the regularizations and constraints we propose in order to enforce prior information on the objects to be reconstructed. Then in section 3, we describe the iterative reconstruction method, *i.e.* the hierarchical optimization procedure based on the closed-form expression of our proximal operator. In section 4, we show that our approach enables quantitative phase estimation of micrometric objects, based on reconstructions of in-line holograms simulated using Lorenz-Mie theory (LMT). Finally, we discuss the quality of the reconstructions based on experimental in-line holograms obtained in two different applications: fluid mechanics with the tracking of an evaporating droplet (size ∼ 100μm) and microscopy of bacteria (size ∼ 1μm).

## 2. Inverse problems formulation

#### 2.1. Hologram formation model

The first step in any inverse problems formulation is defining the image formation model (i.e., the direct model), to describe how given objects produce a hologram on the sensor under incident lighting. Under the hypothesis that the objects of interest are optically thin and located at the same distance from the sensor, they can be modeled by a complex transmittance plane *ṯ*(*x*, *y*) (complex-valued functions or quantities are underlined throughout this article). The amplitude and phase of the transmittance *ṯ*(*x*, *y*) comprise the information to be retrieved from the hologram. The object transmittance plane is illuminated by an incident wave of amplitude $\underset{\_}{{A}_{0}}(x,y)$ and wavelength *λ*. Following the Huygens-Fresnel principle, the object’s plane produces a diffracted wave with a complex amplitude in the hologram plane given by:

*z*is the distance between the sensor and the object plane, $\underset{\_}{{h}_{z}}$ is the propagation convolution kernel and * is the 2-D convolution operator. This kernel is derived from the Rayleigh-Sommerfeld integrale [26] for which it writes:

*i.e.*,

*z*

^{3}≫

*πl*

^{4}/(64

*λ*) where

*l*stands for a characteristic size of the object). The kernel $\underset{\_}{{h}_{z}}$ then writes:

Assuming that the incident beam $\underset{\_}{{A}_{0}}(x,y)=\underset{\_}{A}$ is a plane wave, and hence uniform in the field of view, the intensity of the hologram on the detector is defined by:

The numerical recording of the diffracted light by the camera induces a sampling of the hologram intensity. The acquired data is then defined as a vector ** d** ∈ ℝ

*, with*

^{N}*N*the number of pixels of the sensor. The same discretization can be used in the object plane, the discretized transmittance is then a complex-valued vector

**∈ ℂ**

*ṯ**(the unknowns in our problem). If we represent the 2-D discrete convolution operator by the*

^{N}*N*×

*N*complex-valued matrix

**, the diffraction pattern model**

*H̱**m*at the

*k*-th pixel of the detector is given by:

*indicates the*

_{k}*k*-th element of a vector. This leads to the following data formation model:

*c*that includes the reference plane wave intensity |

**|**

*A̱*^{2}and the conversion factor of the sensor. The vector

**stands for the noise (detection and electronic noise, modeling errors). In the following, the noise component is assumed to be approximately Gaussian and independent on each pixel**

*n**k*.

#### 2.2. Regularized inversion

The inverse problem consists in estimating the transmittance ** ṯ** in the objects plane from a hologram

**. Simple inversion of the hologram formation model Eq. (6) is not possible (**

*d**N*complex-valued unknowns while the hologram has only

*N*real values). It is essential to introduce constraints in addition to requiring that the reconstructed transmittance leads to a hologram model that matches the actual hologram well. These constraints can take two forms: (i) hard constraints that are enforced on all reconstructed transmittances, and (ii) soft constraints that favor solutions with desired properties (e.g., spatial smoothness). Hard constraints restricts feasible transmittances to the set Λ that fulfill the constraints. Soft constraints take the form of a penalization term

*ℛ*. The estimation of the transmittance

**can then be turned into a minimization problem of the form:**

*t**𝒟*corresponds to the opposite of the log-likelihood under our assumption of uncorrelated Gaussian noise:

**is a diagonal matrix with the**

*W**k*-th diagonal entry corresponding to the inverse of the noise variance at pixel

*k*([

**]**

*W*_{k,k}= 0 in the absence of measurement) [27–30]. The scalar factor

*c*is generally not useful, but it is important to set it in order to match the scale of the data

**while keeping a normalized transmittance**

*d***between 0 and 1, as described in section 2.2.2. The minimization problem Eq. (7) takes the classical form of a**

*ṯ**maximum a posteriori*(MAP) estimation, with

*ℛ*the opposite of the log prior and a prior probability ℙ(

**∉ Λ) = 0.**

*ṯ*### 2.2.1. Estimation of the optimal scaling parameter *c*^{*}

The optimal scaling parameter *c*^{*} is given by minimizing the data term Eq. (8) with respect to the unknown *c* :

*c*by

*c*

^{*}(

**) in Eq. (8) defines a data term that depends only on**

*ṯ***:**

*ṯ*### 2.2.2. Physical constraints on the transmittance of the objects

Several physically-grounded constraints on the transmittance of the objects can be considered. Remember that the transmittance ** ṯ** is complex-valued, so let us introduce the notation ℜ𝔢(

**) and ℑ𝔪(**

*ṯ***) which stands for the real and imaginary parts of**

*ṯ***, respectively. Let us also introduce another possible notation for**

*ṯ***as a decomposition in a modulus |**

*ṯ***| and a phase term**

*ṯ**φ*(

**) ∈]−**

*ṯ**π*,

*π*] as follows:

Table 1 lists the constraints we use in this paper. The first two constraints are soft constraints that favor transmittances such that many pixels are fully transparent (|[** ṯ**]

*| = 1 for many*

_{k}*k*) and/or introduce no phase shift (ℑ𝔪([

**]**

*ṯ**) = 0, thus*

_{k}*φ*([

**]**

*ṯ**) = 0 for many*

_{k}*k*). The third constraint favors transmittances that are spatially smooth, while preserving sharp edges (the boundaries of the objects). It uses a generalization of the total-variation to complex-valued images and includes a parameter

*∊*that help prevent a staircasing effect:

*k*-th pixel and the value at the pixel immediately afterward in the

*x*direction, and ${\mathbf{\Delta}}_{k}^{x}$ performs the same kind of difference in the

*y*direction.

The last two constraints are hard constraints, that is to say that they are fully enforced. The first hard constraint prevents the reconstruction of a transmittance with a modulus larger than 1 (no amplification by the medium). Notice that this last constraint, as well as the soft constraint *ℛ*_{1}, requires that ** ṯ** be normalized between 0 and 1 by adjusting the scaling factor

*c*. The second hard constraint can be used to restrict the set of possible phase shifts produced when the incident wave crosses the transmittance plane.

Simple considerations on the sample thickness and maximum contrast between the refractive index of the medium and of the object lead to a bound on the minimum / maximum phase shifts.

The minimization problem including the various physical constraints just discussed is summarized in the next equation:

*α*

_{1},

*α*

_{2},

*α*

_{3}} ∈ ℝ

^{+}are so-called hyper-parameters that balance the importance of each physical constraint.

## 3. Proposed iterative reconstruction algorithm

In the previous section, we gave a formulation of the reconstruction problem as an optimization problem. Solving the optimization Eq. (14) is challenging given the large number of unknowns involved (as many as the number of pixels in the hologram), the non-differentiability of terms *ℛ*_{1} and *ℛ*_{2}, the non-linearity of the constraint $\underset{\_}{\mathit{t}}\in \mathrm{\Omega}\iff \sqrt{\Re \U0001d522{\left(\underset{\_}{\mathit{t}}\right)}^{2}+\Im \U0001d52a{\left(\underset{\_}{\mathit{t}}\right)}^{2}}\le 1$, and the non-convexity of terms *𝒟* and *ℛ*_{1}. To tackle this problem, we use a variable splitting strategy, *i.e.* we introduce additional variables ** a** ∈ ℝ

*and*

^{N}**∈ ℝ**

*b**such that at convergence*

^{N}**=**

*ṯ***+**

*a**i̱*

**. We then recast the optimization problem Eq. (14) as an equivalent constrained optimization problem:**

*b**ℒ*and search for a saddle point:

_{β}**is the scaled dual variable, and**

*u̱**β*> 0 is the augmented Lagrangian penalization parameter, see [31]. To find a saddle point, two steps can be repeated: minimization with respect to the primal variables and update of the dual variables:

*j*corresponds to the iteration number.

Mourya *et al.* [25] showed that optimization Eq. (17) could be performed efficiently by hierarchical optimization when minimizing with respect to some of the variables could be performed in closed form. Following this approach, the algorithm (*cf.* Eq. (17)–(18)) can be turned into two steps: (i) the minimization of a (non-convex) smooth function, and (ii) a dual update:

Evaluations of the cost function in Eq. (19) and its gradient require finding the optimal values {*a*^{*}(** ṯ**),

*b*^{*}(

**)}. Given the non-linear constraints, the non-differentiability and non-convexity of the cost function, this is not trivial. However, the problem is separable and can thus be solved independently at each pixel. In appendix A, we show that a closed-form expression of the global minimum can be obtained. This closed-form expression is the key to the simplicity and efficiency of our iterative reconstruction algorithm.**

*ṯ*The minimization problem Eq. (19) is non-convex but differentiable, a local minimum can then be found using a quasi-Newton algorithm. Here we use the Variable Metric Limited Memory (VMLM) algorithm [32], including a variant for dealing with linear bound constraints (VMLM-B) [33].

In regularized inverse problems approaches, the difficulty of tuning the hyper-parameter for each *a priori* is a central problem which is more complex in this approach that includes three regularization hyper-parameters (*α*_{1},*α*_{2},*α*_{3}), and 2 bound constraints on the phase (*φ*_{min},*φ*_{max}) which are needed to be tuned. Experience shows that the tuning of regularizations requires a tradeoff between the quality of the recovered signal in terms of SNR while keeping a good spatial resolution. As it is complicated to get precise quantitative criterions to evaluate the quality of the reconstruction when the inverse problem uses non-linear regularizations and constraints (*e.g.* Cramer-Rao lower bounds are inapplicable in this context), the quality criterion is most of the time empirical.

Moreover another critical point is the automatic tuning of the hyper-parameters. Different tests showed that a good combination tuned by “hand” for one hologram can be applied to other holograms recorded in the same configuration. This is due to the normalization of the transmittance (its modulus is always in [0, 1]) that make the scale of the spatial gradients and L1 norms comparable from one reconstruction to another. To further simplify the tuning of these hyper-parameters, we normalize each hologram by its maximum value so that the range of the data term is similar from one acquisition to another. There are some special cases that suggest some particular settings of the hyper-parameters:

- to reconstruct purely absorbing objects, a very large value can be used for the hyper-parameter
*α*_{2}(*α*_{2}→ ∞), - to reconstruct purely phase objects, a very large value can be set for the hyper-parameter
*α*_{1}(*α*_{1}→ ∞), - when the difference of refractive index
*n*_{0}−*n*_{object}≤ 0, we can set*φ*_{min}= 0, and if*n*_{0}−*n*_{object}≥ 0 then*φ*_{max}= 0. - to produce a smoother reconstruction of the phase and modulus, the value of the hyper-parameter
*α*_{3}should be increased.

The parameter *∊* is used to prevent from producing staircase artifacts, it can be set between 10^{−4} and 10^{−3}. Using a very small value for *∊* gives piecewise-constant reconstructions - it tends to an exact TV regularization - and significantly slows down the convergence of the algorithm (by up to one order of magnitude).

With the augmented Lagrangian methods, the parameter *β* controls the convergence rate of the algorithm. For convex problems, Mourya *et al.* [25] have shown that the value of *β* has a very limited impact on the convergence speed. Our optimization problem being non-convex, the value of *β* may also impact the local optimum obtained by the algorithm. In practice, we have found that values for *β* in the range [5, 1000] lead to a satisfying convergence speed and very close reconstruction results.

## 4. Results

#### 4.1. Application in fluid mechanics

### 4.1.1. Context and experiments

The reconstruction procedure is tested in the field of fluid mechanics. The experiments consist in imaging droplets of diethyl ether, evaporating in a turbulent flow. The experimental setup comprises a droplet generation device, a turbulence generation and control device and an in-line holographic arrangement, as described in reference [34–36], using the experimental set up Fig. 1(a). The droplets, with a diameter of ∼100 μm, are illuminated by a monochromatic divergent beam of wavelength *λ* = 532 nm. In-line holograms of these droplets are recorded on the CMOS sensor of a Phantom V611 high speed camera at a framerate of 3 kHz. The sensor is composed of 1280 × 800 pixels with a pixel pitch of 20 μm. The distance *z* from the sensor plane to the droplets is about 630 mm. Note that the beam divergence introduces a magnification ratio of about 1.5, which varies weakly with the distance *z*. The holograms of the droplets are then equivalent to holograms recorded with a plane wave illumination for droplet sizes of about 150 μm and a recording distance of about 1 m. Figure 1(b) shows an experimental hologram obtained at one stage of the acquisition.

### 4.1.2. Simulations of evaporating droplets using the Mie theory

In the context of the fluid mechanics experiments described above, our reconstruction method was first tested on synthetic holograms simulated by Lorenz-Mie Theory (LMT), simulating holograms of fluid particles. LMT provides a rigorous solution to the problem of light scattering by a perfect homogeneous and isotropic sphere [37, 38]. This solution has been extended to the case of radially inhomogeneous spheres [39,40]. This extended formulation can be used to simulate the hologram of a liquid droplet surrounded by a vapor cloud [35], at any distance *z* from the droplet. Synthetic images were generated from LMT simulations, taking into account the integration over sensor pixels. The spherical droplet was defined by a radius *r _{d}* = 50 μm and a complex refractive index

*n*= 1.35. The droplet is illuminated by a monochromatic plane wave of wavelength

_{d}*λ*= 532 nm. The vapor cloud, assumed to be spherical, was defined by a refractive index

*n*(

_{v}*r*),

*r*being the radial distance in spherical coordinates. In the following, we simply used an exponential decay from

*n*= 1 + 10

_{s}^{−4}(at the droplet surface) to 1 far from the droplet, according to the formula

*n*(

_{v}*r*) − 1 = (

*n*− 1) exp (−(

_{s}*r*−

*r*)/

_{d}*σ*) where the decay parameter

*σ*= 100 μm. The sensor is composed of 1000 × 1000 pixels with a pixel pitch of 20 μm. The recording distance between the droplet and the sensor is

*z*= 0.5 m. In the vapor cloud, the refractive index is close to 1 and varies slowly. Thus, it can be considered mostly as a phase object. On the contrary, the liquid droplet behaves as an amplitude object (nearly opaque). The complex diffracting object composed of the droplet and the vapor cloud is then a phase-shifting and absorbing object and the reconstructed transmittance is expected to be complex-valued. For the propagation distance considered in this simulation (

*z*= 0.5 m), simulated holograms of a transparent (

*n*= 1.35) and an absorbing sphere (

_{d}*n*= 1.35 − 0.1

_{d}*i̱*) are found to be identical, because the light refracted through the droplet diverges faster and is negligible compared to the diffracted light. Thus, at large distances, both absorbing and transparent droplets can be considered as absorbing objects. Figure 2(a) shows the synthetic hologram, and Fig. 2(b) shows the phase (on a zoomed field of view) of the back-propagation obtained from the fully complex hologram - with the phase - which is considered as the ground truth phase image. Two noisy holograms are generated, adding a uniform and independent Gaussian noise with a signal to noise ratio (SNR) of 100 and 50, respectively. At the location of the central droplet, as it is a purely absorbing object, the values of the phase are inconsistent. Hence we mask these values to 0 to better concentrate on the vapor cloud’s phase. The mask operates on all pixels whose modulus is lower than 0.5.

The reconstructions were performed with our method involving 16 dual updates with 30 iterations between each update to ensure the good convergence of the algorithm. The 5 hyper-paremeters (*α*_{1}, *α*_{2}, *α*_{3}, *∊*, *β*) were tuned “by hand” to obtain a convincing result. The calculation time was about 5 minutes.

For comparison, a reconstruction of the simulated hologram with SNR = 100 was performed with an Error-Reduction Fienup’s algorithm [6] including a support constraint at each iteration. The support is supposed to be known and is the disk of radius 6*σ* = 600 μm (recall that *σ* is the decay parameter of the refractive index radial profile of the simulated vapor cloud). At this distance, the phase shift is less than 1 % of the maximum phase shift at the interface between the droplet and the vapor cloud, so it is negligible because it is hidden in the noise.

Figure 3 compares the reconstruction of the phase with the ground truth phase of Fig. 2(b). We show the reconstructed phases (a,c,e) and the error maps (b,d,f) for the 2 data holograms at SNR = 100 (a,b) and 50 (c,d). The reconstruction at SNR = 100 with Fienup’s algorithm is shown in Fig. 3(e,f). The error maps are defined as the absolute difference between the reconstruction and the ground truth. As already said for the ground truth, we mask the pixels corresponding to the central droplet which is considered as a purely absorbing object, thus leading to inconsistent phase values. Masking operation is done in the same way as for the ground truth image, *i.e.* on all pixels whose modulus is lower than 0.5. We chose not to show the results of reconstructions of the modulus so as to concentrate on analysis of phase quantification.

At first glance, the reconstructions of the phase are close to the ground truth, showing the correct expected range of values. The same observation is made on the reconstructions of the modulus. The error maps allow a finer comparison, and we logically observe that the errors are more important for data SNR = 50. At SNR = 100, our method gives a much better phase estimation than Fienup’s reconstruction. The error map of this reconstruction shown in Fig. 3(f) illustrates the higher errors - the display range of error values has been enlarged to show how far the error increases. This can be due to the fact that in this simulation, the support is too large to eliminate properly the twin-image artifacts involved by the central droplet. Thus Fienup’s algorithm fails in properly separating phase and absorption information. In the case of experimental holograms of evaporating droplets (see section 4.1.3), this problem will be enforced by the fact that the support of the vapor cloud will be very difficult to estimate precisely because we have no a priori information about its shape. For reconstructions of holograms of bacteria (see section 4.2.2), Fienup’s algorithm could perform better because the phase and absorption information is located at the same place. However the problem of estimating correctly the support of unknown-shaped bacteria would remain challenging. On the contrary, our method seems to distinguish properly phase and absorption information even if it is spacially separated. Moreover, sparse regularizations act as an adaptive estimation of the support of the objects both in phase and amplitude.

To better quantify the errors, we calculate 2 metrics EMAX and RMSE that correspond to the normalized maximum error and the normalized root mean squared error, respectively. Normalization is obtained by dividing the absolute maximum error and root mean square error by the maximum value of the ground truth. These quantitative results confirm the previous observations: a lower SNR yields larger errors. EMAX is very important with values around 15%, while RMSE goes from 1.53% for SNR = 100 to 1.68% for SNR = 50, showing that the global phase estimation is very accurate. Quantified errors at SNR = 100 for Fienup’s reconstruction are much higher with EMAX = 49.41% and RMSE = 10.67%, which confirms our previous observations.

### 4.1.3. Reconstructions from the sequence of experimental holograms

Figure 4(a), shows 3 holograms of the same droplet at successive time steps, extracted from an evaporation time sequence (See Visualization 1). The holograms show the typical circular fringes associated to the droplet and the trace of a vapor wake, whose orientation changes with the flow direction. The inversion results are also presented in Fig. 4(b,c) with focus on the reconstructed phase (see Fig. 4(b)) and modulus (see Fig. 4(c)). These results do not quantitatively validate the reconstruction method but illustrate the potential of the method for this particular application. Note however that the order of magnitude of the maximum phase shift, close to the droplet can be estimated by *φ* = *l* * *n _{s}* * 2

*π*/

*λ*, where

*l*is an estimation of the wake depth and

*n*is the refractive index of the ether vapor / air mixture at the droplet surface. Considering

_{s}*l*equals twice the droplet diameter (about 200 μm) and

*n*≈ 1.5 + 10

_{s}^{−4}(based on measurements made on millimeter sized suspended drops [41]) the maximum phase shift is in the order of 0.35 rad. In [36], droplet hologram sequences were reconstructed using an alternative method [27,28]. This parametric method provides highly accurate quantitative measurement of droplet size and position. Note that the

*z*position estimated by this method is taken as the focus distance for calibrating our propagation kernel $\underset{\_}{{h}_{z}}$ (see (2) and (3)). However it is limited to qualitative information on the vapor wake. In the context of fluid mechanics, measuring and tracking both the droplet and the vapor wake is crucial to understand the role played by turbulence in the evaporation of the droplets and more specifically, some deviations observed in [36] in comparison with the evaporation rates predicted by quasi-steady evaporation and drag force modeling. In this context, the results obtained with the proposed method are promising for quantifying both the droplet and the vapor wake.

From a qualitative point of view, Fig. 5 shows that the pseudo hologram data obtained from the repropagation of the reconstruction (see Fig. 5(b)) has a high fidelity with the experimental hologram (see Fig. 5(a)), while reducing a significant proportion of the experimental noise (sensor noise, outliers in the other planes, *etc.*).

Reconstructions in Fig. 4 are performed with our method involving 16 dual updates with 50 iterations between each update to ensure the good convergence of the algorithm. The calculation has been launched on a CPU Intel Core i7-3630QM (2,40 GHz). The calculation time is approximately 15–20 minutes. Figure 6 shows the phase reconstruction of the hologram at the top-left in Fig. 4 at different iteration steps. If the best results are given for the optimization strategy used for the reconstructions shown in Fig. 4 (see Fig. 6(c)), we can see that we have pretty good reconstructions in almost a few minutes after less iterations and/or updates (see Fig. 6(a,b)), even if the convergence is not reached in these cases.

#### 4.2. Application in defocused microscopy for micro-biology

### 4.2.1. Context and experimental bench

The reconstruction was also tested in the field of micro-biology. The experiment consists in imaging stained bacteria with a microscope setup (see Fig. 7). The objects are bacteria that are either gram negative rods (Escherichia coli, noted *E. coli*) or gram positive cocci (Staphylococcus epidermidis, noted *S. epidermidis*). A clonal population of both of these bacteria is isolated on a Petri plate, then both of the isolates are smeared on a microscope slide. The bacteria are then Gram-stained using standard procedure (bioMérieux PreviColor Gram system). The resulting objects immobilized on the surface of the slide are pink (gram negative) or violet (gram positive), with a typical size of ∼ 1*μm*.

These objects are then imaged using two different microscopic setups (see Fig. 7(a)):

- The first (“reference”) setup is similar to the one typically used in routine in-vitro diagnotics laboratories. It consists in a white-light, transmission, oil-immersion microscope setup, with a color camera, mounted in a standard instrument (Olympus BX-61). The light source is a white LED (5500 K, Mightex FCS-0000-000), coupled with a 600
*μm*-core optical fiber bundle. The end of the optical fiber is positioned in a distance of a few centimeters from the slide, which is placed on a motorized stage. The objects are positioned by the operator at the focal plane of the objective. The image is collected through a 60X/NA1.4 objective (Olympus), a 10× tube lens, and a 5Mpix color camera (Basler, pixel size 3.45*μm*). - Using the same mount, a second setup can be implemented. Using another fiber of the bundle, the light source is switched to a monochromatic LED source (617
*nm*, Mightex FCS-0617-000). Using the motorized stage, the objects are positioned at 32 μm distance from the focal plane of the objective. A band-pass filter (610*nm*, 5*nm*width) is positioned between the objective and the tube lens. This allows acquiring an in-line hologram of the objects (see Fig. 7(b)). For sake of comparisons, images of the object positioned in the focal plane are also acquired with this mount (see Fig. 8).

### 4.2.2. Reconstructions from holograms of bacteria

Figure 7(b) shows a defocused hologram at *z* = 32 μm of a biological sample containing the two species of bacteria, *S. epidermidis* and *E. coli* (*E.coli*). Figure 8(a) shows the white-light illuminated intensity image of the biological sample at focus which stands for the reference image. Figure 8(b) shows the intensity image, at the focal plane of the objective, of the biological sample at the illumination wavelength *λ* = 610 nm. Figure 9(a,b) respectively show the modulus and phase of the reconstruction obtained with our method from the defocused hologram of Fig. 7(b). In the experiment, bacteria are Gram stained, in pink color (chromophore Safranin) for *E.coli* bacteria (gram-negative), and violet color (chromophore Crystal Violet) for *S. epidermidis* (gram-positive). For the 610 nm illumination wavelength, *E.coli* does not absorb and behaves mostly as a phase object, while *s. epidermidis* both absorbs light and shifts the phase. These results are coherent with the expected absorptions of the chromophores for the wavelength 610 nm (see Fig. 11).

Acquiring defocused intensity images (holograms) brings information from phase objects that can not be recorded at the focus plane. From the field of view shown in Fig. 9, we extract 3 regions of interest A, B, and C. Figure 10 show zooms on these 3 regions, for each reference and reconstruction images (see Figs. 8 and 9). We can observe that *E.coli* bacteria are reconstructed by our method like a purely phase object (third and fourth columns), whereas *E.coli* bacteria cannot be observed on the “red”-light intensity image at focus (second column). *S. epidermidis* bacteria are reconstructed by the proposed algorithm like phase and absorbing objects (third and fourth columns). These results show that phase information is promising to discriminate the two populations involved in this experiment, because their phase and modulus characteristics are different. Hence our reconstruction algorithm demonstrates its ability to correctly retrieve phase and absorption information, in accordance with what is physically expected from Fig. 11: a bacteria colored by crystal violet will absorb red light at the wavelength 610 nm while another that is colored by safranin will absorb very little light. No contrast is expected in terms of absorption. Only the phase shift produced by the bacteria makes it visible in a defocused image. This fits with our reconstruction where bacteria colored by crystal violet are visible in the reconstructed modulus, while safranin-colored bacteria are visible only in the reconstructed phase.

## 5. Conclusion & discussion

We have devised a regularized inverse problems approach to quantitatively reconstruct the transmittance of absorbing and phase objects from a single in-line hologram. The originality of the proposed method is the inclusion of several physical constraints and the derivation of an efficient algorithm based on the closed-form expression of the global minimum of a non-convex, non-smooth optimization sub-problem that involves non-linear constraints. From a robustness point of view, this approach was successfully applied on several in-line holography configurations for different applications: fluid mechanics and micro-biology. The performance of the phase reconstruction has been demonstrated both quantitatively with simulated holograms and qualitatively with experimental holograms for applications in fluids mechanics and micro-biology. Particularly, our method have demonstrated its ability to correctly separate phase and absorption information in a coherent way, showing at the same time that phase information is of primary importance for the targeted applications.

A critical point of our method is the large calculation time. For the moment, our algorithm has involved no optimization for its implementation, as it has been initially developed on a single CPU. Other works in this field have shown that GPU acceleration can drastically speed up reconstructions, using fast GPU implementations of modelization and optimization tools classically used in digital holography [42], or applying GPU-accelerated deep neural networks [12]. An implementation on GPU is also possible for our algorithm, as many parts of it can be highly parallelized: (i) the forward model which uses Fast Fourier Transforms ; (ii) the proximal operator which is separable on each pixel of the image.

A limitation could be that the reconstruction is based on a hologram formation model based on the 2-D projection of 3-D objects, which is only valid for thin objects. The proposed inverse problem approach can be easily adapted to reconstruct several absorbing and phase objects located in different 2-D planes at various distances *z* from the sensor plane.

## A. Appendix A

Closed-form expression of the solution to problem (20) pixel by pixel (*s* and *b* are scalar):

The proximal operator

- in the case where the constraints are not active, it corresponds to the solution of the following problem:$$\underset{\mathit{a},\mathit{b}}{\text{arg min}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\beta {\left(p-a\right)}^{2}+\beta {\left(q-b\right)}^{2}+{g}_{{\alpha}_{1},{\alpha}_{2}}(a,b)$$which writes:$$\left({a}^{*},{b}^{*}\right)=\left(p+\frac{{\alpha}_{1}}{2\beta}\frac{p}{\left|p+\underset{\_}{i}{S}_{{\alpha}_{2}/2\beta}(q)\right|},{S}_{{\alpha}_{2}/2\beta}(q)+\frac{{\alpha}_{1}}{2\beta}\frac{{S}_{{\alpha}_{2}/2\beta}(q)}{\left|p+\underset{\_}{i}{S}_{{\alpha}_{2}/2\beta}(q)\right|}\right)$$where
*S*_{α2/2β}(*q*) corresponds to the soft-thresholding of*q*:$${S}_{{\alpha}_{2}/2\beta}(q)\equiv \{\begin{array}{lll}q-\frac{{\alpha}_{2}}{2\beta}\hfill & \text{if}\hfill & q\ge \frac{{\alpha}_{2}}{2\beta}\hfill \\ 0\hfill & \text{if}\hfill & q\in \left[-\frac{{\alpha}_{2}}{2\beta},\frac{{\alpha}_{2}}{2\beta}\right]\hfill \\ q+\frac{{\alpha}_{2}}{2\beta}\hfill & \text{if}\hfill & q\le -\frac{{\alpha}_{2}}{2\beta}\hfill \end{array}$$Special case: if p=0 and q=0 - if the constraint on the phase is active, then
*φ*(*a*^{*}+*i̱b*^{*}) =*φ*or_{min}*φ*(*a*^{*}+*i̱b*^{*}) =*φ*and the associated modulus |_{max}*a*^{*}+*i̱b*^{*}| gives:$$\left|{a}^{*}+\underset{\_}{i}{b}^{*}\right|=\{\begin{array}{ll}1\hfill & \text{if}\hspace{0.17em}{\rho}^{*}({\phi}^{*})\ge 1\hfill \\ {\rho}^{*}({\phi}^{*})\hfill & \text{if}\hspace{0.17em}\hspace{0.17em}{\rho}^{*}({\phi}^{*})\in [0,1]\hfill \\ 0\hfill & \text{if}\hspace{0.17em}{\rho}^{*}({\phi}^{*})\le 0\hfill \end{array}$$with$${\rho}^{*}(\phi )=p\text{cos}(\phi )+q\text{sin}(\phi )+\frac{{\alpha}_{1}}{2\beta}-\frac{{\alpha}_{2}}{2\beta}\left|\text{sin}(\phi )\right|$$and$$\begin{array}{ll}{\phi}^{*}=\hfill & \underset{\phi \in \{{\phi}_{\mathit{min}},{\phi}_{\mathit{max}}\}}{\text{arg min}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\beta {\left(p-{\rho}^{*}(\phi )\text{cos}(\phi )\right)}^{2}\hfill \\ \hfill & \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+\beta {\left(q-{\rho}^{*}(\phi )\text{sin}(\phi )\right)}^{2}\hfill \\ \hfill & \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\alpha}_{1}\left(1-{\rho}^{*}(\phi )\right)\hfill \\ \hfill & \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\alpha}_{2}{\rho}^{*}(\phi )\left|\text{sin}(\phi )\right|.\hfill \end{array}$$

## Funding

Centre National de la Recherche Scientifique; La Région Auvergne-Rhône-Alpes; Agence Nationale de la Recherche (ANR-11-LABX-0063, ANR-11-IDEX-0007).

## Acknowledgments

This work has been supported in part by the CNRS under project RESSOURCES (Défi Imag’In). Funding for this project was also provided by a grant from “La Région Auvergne-Rhône-Alpes”. This work was performed within the framework of the LABEX PRIMES (ANR-11-LABX-0063) of Université de Lyon, within the program “Investissements d’Avenir” (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR).

## References and links

**1. **D. Gabor, “A new microscopic principle,” Nature **161**, 777–778 (1948). [CrossRef] [PubMed]

**2. **M. Liebling, T. Blu, and M. Unser, “Complex-wave retrieval from a single off-axis hologram,” J. Opt. Soc. Am. A **21**, 367–377 (2004). [CrossRef]

**3. **I. Yamaguchi, J.-i. Kato, S. Ohta, and J. Mizuno, “Image formation in phase-shifting digital holography and applications to microscopy,” Appl. Opt. **40**, 6177–6186 (2001). [CrossRef]

**4. **M. Jericho, H. Kreuzer, M. Kanka, and R. Riesenberg, “Quantitative phase and refractive index measurements with point-source digital in-line holographic microscopy,” Appl. Opt. **51**, 1503–1515 (2012). [CrossRef] [PubMed]

**5. **R. W. Gerchberg and W. O. Saxton, “A practical algorithm for the determination of the phase from image and diffraction plane pictures,” Optik **35**, 237 (1972).

**6. **J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. **21**, 2758–2769 (1982). [CrossRef] [PubMed]

**7. **R. Horisaki, Y. Ogura, M. Aino, and J. Tanida, “Single-shot phase imaging with a coded aperture,” Opt. Lett. **39**, 6466–6469 (2014). [CrossRef] [PubMed]

**8. **Z. Wang, Q. Dai, D. Ryu, K. He, R. Horstmeyer, A. Katsaggelos, and O.S. Cossairt, “Dictionary-based phase retrieval for space-time super resolution using lens-free on-chip holographic video,” in *Imaging and Applied Optics 2017 (3D, AIO, COSI, IS, MATH, pcAOP)*, OSA Technical Digest (online) (Optical Society of America, 2017), paper CTu2B.3. [CrossRef]

**9. **Y. Rivenson, Y. Wu, H. Wang, Y. Zhang, A. Feizi, and A. Ozcan, “Sparsity-based multi-height phase recovery in holographic microscopy,” Sci. Rep. **6**, 37862 (2016). [CrossRef] [PubMed]

**10. **D. Ryu, Z. Wang, K. He, G. Zheng, R. Horstmeyer, and O. Cossairt, “Subsampled phase retrieval for temporal resolution enhancement in lensless on-chip holographic video,” Biomed. Opt. Express **8**, 1981–1995 (2017). [CrossRef] [PubMed]

**11. **F. Eilenberger, S. Minardi, D. Pliakis, and T. Pertsch, “Digital holography from shadowgraphic phase estimates,” Opt. Lett. **37**, 509–511 (2012). [CrossRef] [PubMed]

**12. **Y. Rivenson, Y. Zhang, H. Günaydın, D. Teng, and A. Ozcan, “Phase recovery and holographic image reconstruction using deep learning in neural networks,” arXiv preprint arXiv:1705.04286 (2018).

**13. **A. Sinha, J. Lee, S. Li, and G. Barbastathis, “Lensless computational imaging through deep learning,” Optica **4**, 1117–1125 (2017). [CrossRef]

**14. **S. Sotthivirat and J. A. Fessler, “Penalized-likelihood image reconstruction for digital holography,” J. Opt. Soc. Am. A **21**, 737–750 (2004). [CrossRef]

**15. **L. Denis, D. Lorenz, E. Thiébaut, C. Fournier, and D. Trede, “Inline hologram reconstruction with sparsity constraints,” Opt. Lett. **34**, 3475–3477 (2009). [CrossRef] [PubMed]

**16. **A. Bourquard, N. Pavillon, E. Bostan, C. Depeursinge, and M. Unser, “A practical inverse-problem approach to digital holographic reconstruction,” Opt. Express **21**, 3417–3433 (2013). [CrossRef] [PubMed]

**17. **F. Soulez, É. Thiébaut, A. Schutz, A. Ferrari, F. Courbin, and M. Unser, “Proximity operators for phase retrieval,” Appl. Opt. **55**, 7412–7421 (2016). [CrossRef] [PubMed]

**18. **C. Schretter, D. Blinder, S. Bettens, H. Ottevaere, and P. Schelkens, “Regularized non-convex image reconstruction in digital holographic microscopy,” Opt. Express **25**, 16491–16508 (2017). [CrossRef] [PubMed]

**19. **J. Song, C. L. Swisher, H. Im, S. Jeong, D. Pathania, Y. Iwamoto, M. Pivovarov, R. Weissleder, and H. Lee, “Sparsity-based pixel super resolution for lens-free digital in-line holography,” Sci. Rep. **6**, 24681(2016).

**20. **S. Bettens, H. Yan, D. Blinder, H. Ottevaere, C. Schretter, and P. Schelkens, “Studies on the sparsifying operator in compressive digital holography,” Opt. Express **25**, 18656–18676 (2017). [CrossRef] [PubMed]

**21. **A. Berdeu, F. Momey, B. Laperrousaz, T. Bordy, X. Gidrol, J.-M. Dinten, N. Picollet-D’hahan, and C. Allier, “Comparative study of fully three-dimensional reconstruction algorithms for lens-free microscopy,” Appl. Opt. **56**, 3939–3951 (2017). [CrossRef] [PubMed]

**22. **N. Parikh and S. Boyd, Proximal Algorithms, Foundations and Trends® in Optimization **1**, 127–239 (2014). [CrossRef]

**23. **P. L. Combettes and J.-C. Pesquet, “Proximal Splitting Methods in Signal Processing,” in *Fixed-Point Algorithms for Inverse Problems in Science and Engineering,* (Springer, 2011), pp. 185–212. [CrossRef]

**24. **W. Hare and C. Sagastizábal, “Computing proximal points of nonconvex functions,” Math. Program. **116**, 221–258 (2009). [CrossRef]

**25. **R. Mourya, L. Denis, J.-M. Becker, and E. Thiebaut, “Augmented lagrangian without alternating directions: Practical algorithms for inverse problems in imaging,” in *Proc. of 2015 IEEE International Conference on Image Processing (ICIP)*, (IEEE, 2015), pp. 1205–1209. [CrossRef]

**26. **J. W. Goodman, *Introduction to Fourier Optics*, (Roberts and Company Publishers, 2005).

**27. **F. Soulez, L. Denis, C. Fournier, É. Thiébaut, and C. Goepfert, “Inverse-problem approach for particle digital holography: accurate location based on local optimization,” J. Opt. Soc. Am. A **24**, 1164–1171 (2007). [CrossRef]

**28. **F. Soulez, L. Denis, E. Thiébaut, C. Fournier, and C. Goepfert, “Inverse problem approach in particle digital holography: out-of-field particle detection made possible,” J. Opt. Soc. Am. A **24**, 3708–3716 (2007). [CrossRef]

**29. **C. Fournier, L. Denis, M. Seifi, and T. Fournel, “Digital hologram processing in on-axis holography,” Multi-Dimensional Imaging **0**, 51–73 (2014). [CrossRef]

**30. **C. Fournier, F. Jolivet, L. Denis, N. Verrier, E. Thiebaut, C. Allier, and T. Fournel, “Pixel super-resolution in digital holography by regularized reconstruction,” Appl. Opt. **56**, 69–77 (2017). [CrossRef]

**31. **S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Foundations and Trends® in Machine Learning **3**, 1–122 (2011). [CrossRef]

**32. **J. Nocedal, “Updating quasi-Newton matrices with limited storage,” Math. Comput. **35**, 773–782 (1980). [CrossRef]

**33. **E. Thiébaut, “Optimization issues in blind deconvolution algorithms,” Proc. SPIE **4847**, pp. 174–183 (2002). [CrossRef]

**34. **D. Chareyron, J.-L. Marié, C. Fournier, J. Gire, N. Grosjean, L. Denis, M. Lance, and L. Méès, “Testing an in-line digital holography ‘inverse method’ for the Lagrangian tracking of evaporating droplets in homogeneous nearly isotropic turbulence,” New J. Phys. **14**, 043039 (2012). [CrossRef]

**35. **L. Méès, N. Grosjean, D. Chareyron, J.-L. Marié, M. Seifi, and C. Fournier, “Evaporating droplet hologram simulation for digital in-line holography setup with divergent beam,” J. Opt. Soc. Am. A **30**, 2021–2028 (2013). [CrossRef]

**36. **J.-L. Marié, T. Tronchin, N. Grosjean, L. Méès, O. C. Öztürk, C. Fournier, B. Barbier, and M. Lance, “Digital holographic measurement of the lagrangian evaporation rate of droplets dispersing in a homogeneous isotropic turbulence,” Exp. Fluids **58**, 11 (2017). [CrossRef]

**37. **C. F. Bohren and D. R. Huffman, *Absorption and Scattering of Light by Small Particles*, (Wiley, 1983).

**38. **G. Gouesbet and G. Gréhan, *Generalized Lorenz-Mie Theories*, (Springer, 2017). [CrossRef]

**39. **L. Kai and P. Massoli, “Scattering of electromagnetic-plane waves by radially inhomogeneous spheres: a finely stratified sphere model,” Appl. Opt. **33**, 501–511 (1994). [CrossRef] [PubMed]

**40. **F. Onofri, G. Gréhan, and G. Gouesbet, “Electromagnetic scattering from a multilayered sphere located in an arbitrary beam,” Appl. Opt. **34**, 7113–7124 (1995). [CrossRef] [PubMed]

**41. **G. Toker and J. Stricker, “Holographic study of suspended vaporizing volatile liquid droplets in still air,” Int. J. Heat Mass Tran. **39**, 3475–3482 (1996). . [CrossRef]

**42. **Y. Endo, T. Shimobaba, T. Kakue, and T. Ito, “GPU-accelerated compressive holography,” Opt. Express **24**, 8437–8445 (2016). [CrossRef] [PubMed]