Reconstruction of in-line holograms: combining model-based and regularized inversion.

In-line digital holography is a simple yet powerful tool to image absorbing and/or phase objects. Nevertheless, the loss of the phase of the complex wavefront on the sensor can be critical in the reconstruction process. The simplicity of the setup must thus be counterbalanced by dedicated reconstruction algorithms, such as inverse approaches, in order to retrieve the object from its hologram. In the case of simple objects for which the diffraction pattern produced in the hologram plane can be modeled using few parameters, a model fitting algorithm is very effective. However, such an approach fails to reconstruct objects with more complex shapes, and an image reconstruction technique is then needed. The improved flexibility of these methods comes at the cost of a possible loss of reconstruction accuracy. In this work, we combine the two approaches (model fitting and regularized reconstruction) to benefit from their respective advantages. The sample to be reconstructed is modeled as the sum of simple parameterized objects and a complex-valued pixelated transmittance plane. These two components jointly scatter the incident illumination, and the resulting interferences contribute to the intensity on the sensor. The proposed hologram reconstruction algorithm is based on alternating a model fitting step and a regularized inversion step. We apply this algorithm in the context of fluid mechanics, where holograms of evaporating droplets are analyzed. In these holograms, the high contrast fringes produced by each droplet tend to mask the diffraction pattern produced by the surrounding vapor wake. With our method, the droplet and the vapor wake can be jointly reconstructed.


Introduction
The in-line holography principle was first introduced by Gabor in the late 1940s [1]. As schematically presented in Fig. 1(a), in this kind of microscopy, a coherent light source is scattered by the sample, which produces a diffracted wave. The incident beam and the diffracted wave interfere on the sensor plane. The sensor records the resulting intensity. These interferences are sensitive to the amplitude changes and phase shifts induced by the objects, making this technique particularly suitable for image and study absorbing and/or phase samples. Today, this high sensitivity is used in numerous fields including biology [2][3][4][5][6][7], fluid mechanics [8][9][10][11], and particle characterization [12][13][14].
The simplicity of the experimental setup is counterbalanced by an inherent drawback: the lack of a reference arm. This leads to a loss of the information concerning the wavefront phase on the sensor plane. Based on the solution of the diffraction equations [15], simple backpropagation can be applied to refocus a hologram. Nonetheless, this technique leads to strong artifacts such as border effects and twin-images. Thus, dedicated algorithms need to be developed to improve the reconstruction of the objects [16].
Numerous iterative algorithms proposed in the literature are based on a phase retrieval formulation, where the unknown is the unmeasured phase on the sensor plane. After the phase recovery, the complex amplitude at the sensor plane is backpropagated, leading to improved focusing in the object space [17,18]. An inherent limitation of such approaches is the presence of artifacts due to signal truncation at the borders of the hologram. These approaches also focus on reconstructing the complex wavefront and not the object itself, leading to a non straightforward analysis of the object's properties. Recent techniques have emerged, based on inverse problems formulations, where the unknown is directly the object. The reconstruction consists in maximizing the fidelity between the actual measurements and the data simulated by a hologram formation model, the so-called "forward model".
Two classes of inversion methods exist depending on the object to be reconstructed: • If the object is sufficiently simple, it can be described by a limited number of parameters, such as its 3D location, shape parameters or optical parameters. In that case, the forward model depends only on these parameters, and model fitting approaches, socalled "parametric" approaches, provide robust and highly accurate reconstructions of the samples [6, 9, 12-14, 19, 20]. Model fitting and greedy algorithms belong to this class of methods.
• If the object is too complex to be described by a few parameters, it is approximated in a pixelated space such as a transmittance plane sampled on a 2D grid [21] or a scattering potential sampled on a 3D grid [22]. In that case, the forward model consists in numerically propagating the diffracted wave through the sample, up to the sensor plane. The addition of physical constraints and wisely chosen regularizations is then necessary to solve the illposed problem, in so-called "non-parametric" approaches. Recent works have demonstrated the feasibility of such regularized reconstructions to retrieve the lost phase information from a single in-line hologram of absorbing and/or phase objects [4,23,24]. Penalized-likelihood image reconstruction [25], Maximum A Posteriori [26] and Compressive sensing [27][28][29] frameworks belong to this second class of methods.
In this work, we combine the two approaches by reconstructing the sample using both models, as presented in Fig. 1(b). The aim of combining the two complementary methods is to obtain a more accurate model of the object, leading to better reconstructions.

Proposed method
The objective is to reconstruct an object from its experimental in-line intensity hologram I d , assuming that the object can be divided into two subparts u and v. A diffraction model is used for each of these two parts to predict the diffracted wave U di f on the sensor: a parametric model for u and a non-parametric model for v as seen on Fig. 1(b).
The idea is to combine an accurate parametric reconstruction of simple shaped objects by model fitting and a regularized reconstruction with penalties adapted to the remaining part of the sample. For instance, in the chosen example in fluid mechanics, this approach enables the simultaneous accurate reconstruction of droplets that act as amplitude objects and a vapor wake that acts as a pure phase object. By accounting for two diffraction models, the proposed joint approach aims at improving the reconstruction of all the objects. The global forward model consists in simulating the intensity I s in the hologram plane. It is obtained from the linearity of the wave equation by the superposition of the incident wave U inc and the waves diffracted by the two parts of the object described by u and v: Fig. 1. (a) General scheme of in-line holography: a coherent light source produces a coherent incident wavefront U inc that is scattered by the sample, producing a diffracted wavefront U di f . These wavefronts interfere on the sensor plane that records the resulting Note that optional beam shaping optics or imaging optics to magnify the image are not depicted. (b) Scheme of the proposed model in the case of an evaporating droplet: the object to be reconstructed is composed of two subparts. The wave diffracted by the spherical droplet is given by the Mie model, which is an analytic parametric solution of the diffraction equations. The vapor wake is described by a phase transmittance plane whose diffracted wave is propagated according to the Rayleigh-Sommerfeld theory. Based on the linearity of the equations of the diffraction, these two subparts interfere to create the total wave diffracted by the object: The underlying hypothesis corresponds to the first Born approximation [30]: each part of the sample diffracts an incident wavefront U inc that is assumed to be unperturbed.
Reconstructing u and v can be stated in the form of a problem consisting in minimizing the residues between the data I d and the intensity I s simulated by the model, sampled at the same spatial points as the sensor pixels (i, j). To do so, the estimationsû andv of u and v are obtained by solving the following optimization problem: with the cost function C dat a defined by: where c is a real-valued coefficient to account for the scaling of the acquisition dynamics and W is a matrix of weighting coefficients to possibly account for the presence of defective pixels or the limited size of the sensor [13,31]. These coefficients are inversely proportional to the pixel variances and account for the confidence one has in the data value [I d ] i, j at pixel with 2D coordinates (i, j). For example, for a defective pixel, [W ] i, j = 0 (see appendix for more details).
To solve this joint problem (2), the implemented reconstruction scheme alternates between optimization of the parametric model depending on u and optimization of the non-parametric model depending on v.
The optimization on u consists in estimating a limited number of parameters from a large number of measurements (typically a picture or a sequence of pictures with hundreds of thousands of pixels). As a consequence, the dedicated methods are robust and accurate [12,13].
Conversely, minimizing (2) on v is not sufficient to obtain a satisfactory reconstruction since the image-based problem is ill-posed: several solutions v can solve the minimization problem, by perfectly matching the noisy data. Yet, none of these correspond to a good reconstruction because of the noise amplification phenomenon. It is thus necessary to stabilize the estimation by adding adequate prior information about the object. These priors take the form of hard constraints on the domain D (v) admissible for v and/or regularization terms C r reg in the cost function, whose aim is to favor more realistic solutions. The optimization with respect to v consequently becomes, for a given u: where µ r are hyperparameters to balance the regularization terms C r reg compared to the data fidelity term C dat a .
The proposed method is summarized by algorithm 1.

Algorithm 1 Hologram reconstruction by a joint diffraction model
Optimization with u fixed 7: return (u, v)

Results
We illustrate the reconstruction method with an application in fluid mechanics, on holograms of evaporating droplets surrounded by a vapor plume. It was recently proposed on similar types of datasets to combine a parametric estimation of the spherical objects with a simple backpropagation of the residues [9]. This method succeeds in reconstructing a rough orientation and shape of the vapor plume but fails to provide any information about the object phase. Indeed, as it works directly on differences of intensities, it can not correctly account for the interferences between the different parts of the sample.
With our proposed method, we propose a more accurate forward model that accounts for these interferences, with the aim of obtaining better reconstructions.
In this case, the parametric part of the object is the set of the N d spherical droplets present in the sample. These droplets can be described by a limited number of parameters: respectively their 2D position on the sensor plane, their distance to the sensor and their radius. Due to the large distances between the droplets and the sensor, their refractive index n d has a negligible influence on their hologram and is consequently not a fitted parameter.
The non-parametric part v of the object corresponds to the map of the phase shift introduced by the surrounding medium and the vapor wakes generated by evaporation of the droplets. Under the hypothesis of an optically thin medium, it represents the 2D projection of the 3D object along the line of sight. U is obtained with a multi-z Rayleigh-Sommerfeld propagation (see appendix for more details on the forward model).
A simulation and an experimental hologram are reconstructed. In both cases, N it = 20 iterations are performed in algorithm 1.

Data simulation
In this section, a synthetic hologram with theoretical parameters u th to test the proposed method as presented in Fig. 2(a). The sample is composed of N d = 2 droplets placed at x th 1 = y th 1 = −1.7 mm and x th 2 = y th 2 = 1.7 mm at distances of z th 1 = 0.5 m and z th 2 = 0.45 m. They have a radius of r th 1 = 50 µm and r th 2 = 70 µm. In the simulation, the first droplet is surrounded by a spherical vapor cloud whose refractive index n cloud (r) in terms of radial distance r from the center of the droplet presents an exponential decay from the droplet surface: where n s = n air + 10 −4 is the refractive index at the surface of the droplet and the decay parameter is σ = 100 µm. The cloud is numerically cut at r c = 6σ. The phase shift produced by this vapor cloud is obtained by integration along the line of sight as presented in Fig. 2(b). This projection serves as the ground truth for the non-parametric reconstruction v th . The second droplet is not surrounded by a vapor cloud and consequently does not induce any phase shift in its vicinity.
The hologram of this synthetic object must be simulated. As mentioned in appendix, the Lorenz-Mie theory provides a rigorous solution of the light propagation equations for homogeneous and isotropic spheres. This model has been extended to radially inhomogeneous spheres [32,33]. This makes possible to simulate the hologram produced by an evaporating homogeneous droplet surrounded by a vapor cloud with spherical symmetry [19].
In this simulation, the sensor plane is composed of 512 × 512 pixels with a pitch of p = 20 µm. The illumination is a monochromatic plane wave of wavelength λ = 532 nm with a normalized amplitude of U inc = 1 on the sensor plane.
Finally, a white Gaussian noise of standard deviation σ = 0.02 (SNR = 50) is added to the simulated hologram I th and a scaling factor c th = 0.9 is applied: The results are shown in Fig. 2 and Table 1.

Reconstructed parametric part u
As expected, Figs. 2(a)-2(c) show that the hologram produced only by the parametric part u of the model perfectly matches the signal produced by the second droplet, but fails to explain the bright center of the hologram of the first droplet. This characteristic feature is due to the surrounding vapor cloud [19].
Focusing on the blue values in Table 1, it will be noted that the parameters fitted by the parametric part of the model correctly retrieve the theoretical values. The errors on the x ypositions and the radii are compared to the pixel pitch p = 20 µm. They are less than 5 %, that is to say under the micron in a field of several millimeters and for sizes of several tens of microns. Concerning the propagation distances, the errors oscillate around a few hundreds of microns for theoretical values of several tens of centimeters.
The high precision of the droplet parameters supports the proposed joint method. Indeed, a reconstruction only based on a non-parametric model would have a resolution limited by the pixel size, and consequently an accuracy on the radii and the x y-positioning in the same order of magnitude. In the present case, the parametric model achieves a precision two orders of magnitude below the pixel size. Table 1 shows that the accuracy on the second droplet d = 2 parameters estimation is one order of magnitude better than for the first droplet d = 1. This result was expected as its hologram is not perturbed by a vapor cloud. In addition, the first (in red) and final (in blue) estimates in Table 1 as well as the curves in Fig. 2(f), show that the accuracy of the parametric estimation is not affected by the fact that the algorithm has the degree of freedom to reconstruct a phase in the

Errors on parameters estimation
The fitting of the parameters of the first droplet, whose diffraction pattern is impacted by the vapor cloud, is slightly improved by the joint modeling (see column d = 1 in table 1).

Reconstructed non-parametric part v
Figure 2(d) shows the phase shift reconstructed by non-parametric part of the model. The structure, the shape and the support of the vapor cloud of the first droplet is in good agreement with its theoretical projection as shown in Fig. 2(b). These results support the validity of the optically thin media hypothesis.
Interestingly, the reconstructed vapor cloud is not continuous and is reset to zero on the droplet support by the reconstruction algorithm. This supports the fact that the phase shift information just before and after the droplet is lost through the refraction process and can not be reconstructed.
Given the difference between the reconstruction and the theory, as shown in the insets in Fig. 2(c), it appears that the phase is slightly underestimated, apparently due to the regularizations. Nonetheless, these terms can not be further diminished without producing noisy reconstructions. Thus, the missing signal is under the noise level.
As expected, no phase is reconstructed around the second droplet. This shows that the proposed joint model correctly splits the phase shift effects from the droplets signal. Since the fringes of the second droplets are fitted by the parametric model, no phase is needed in the non-parametric part of the object.
Combining the two parts of the model allows a perfect fit of the simulated hologram as no structure is observed in the residues: only the noise remains in Fig. 2(e).
The chosen regularizations aslo prevent the appearance of any twin-image artifacts around the droplets and the vapor cloud, which are generally present in in-line holographic reconstructions.

Data acquisition
The proposed method is applied to the experimental holograms of diethyl ether droplets evaporating in a turbulent flow. As further described in [9,19,24], a piezoelectric jetting device injects droplets of ∼ 100 µm in a turbulent flow generated and controlled by three pairs of woofers. The in-line holographic setup is composed of a monochromatic beam of wavelength λ = 532 nm and a Phantom V611 high speed camera with a framerate of 3 kHz. The CMOS sensor is composed of 1280 × 800 pixels with a pitch of 20 µm. Figure 3 shows the reconstruction of a hologram extracted from a sequence of 208 holograms acquired over 70 ms (see Visualization 1). In Fig. 3(d), the experimental hologram presents the typical circular fringes with high spatial frequencies associated with the droplets. The smooth traces of three vapor wakes can easily be identified. In addition, a global structured faint signal is present all over the field of view.

Reconstructed parametric part u
The parametric part of the model actually fits eight droplets, one being slightly out of the field of view as sketched in Fig. 3(b). As previously observed in the simulated case, the corresponding intensity, given in Fig. 3(e) only accounts for the high contrast fringes but can not explain the overall background evolution neither the vapor wakes hologram. Fig. 3(b) shows that, despite the fact that their holograms fill extended parts of the field of view, the droplets are so small that they only spread over a few pixels in a non-parametric representation and are pixelated. As previously concluded for the simulated case, this highlights the fact that a parametric approach not based on pixels reconstruction is particularly adapted to explain this part of the sample.
The parameter evolutions along the iterations resembles those emphasized in the simulation analysis of Fig. 2(f) and the parameter refinements also stop evolving after the 10 th iteration. The corresponding curves are consequently not presented here. Table 2 sums up their initial and final values. The diversity of propagation distances z d that are present in the field of view, ranging from 0.257 m to 0.961 m, confirms the relevance of the multi-z approach proposed in the forward model (see appendix). Like the simulated case, the change in the parameters along the iteration is several orders of magnitude below their initial estimates. This can be explained by the fact that the vapor wake perturbs only a small part of the circular fringes. Most of the signal on which the parametric approach bases its reconstruction is consequently unperturbed and the reconstruction accuracy only slightly benefits from the vapor wake estimation. Nonetheless, like the simulated case, the first three droplets that produce a strong vapor wake, present the strongest relative change in their parameter estimates, supporting the fact that even if it is small, there is an improvement in these parameters along the iterations of the alternate approach.

Reconstructed non-parametric part v
The non-parametric part of the reconstructed object, given in Fig. 3(c), clearly shows the three vapor wakes as well as the faint turbulent background. As shown on Fig. 3(f), the corresponding intensity explains part of the high frequencies of the recorded hologram, generated by the sharp edges of the droplets shadow, as well as the large scale faint holograms beyond the reach of the parametric model.
Thus, the vapor wakes spread over long distances while the turbulent background presents large-scale structures. These properties support the choice of a non-parametric approach to reconstruct them. The global forward model accounting for the interferences between the two subparts, it is possible to retrieve the phase information lost in the sensor plane during the acquisition, as presented in Fig. 3(g), by combining the phase predicted by the parametric model, given in Fig 3(h), and the phase propagated by the non-parametric approach, given in Fig. 3(i).
Finally, one may wonder if the background structures reconstructed by the non-parametric approach are artifacts due to the regularization or are an actual signal. In the movie in Visualization 1, the 208 frames of the sequence have been reconstructed independently. The large scale structures appear to evolve coherently frame after frame. They are consequently a meaningful signal present in the data rather than reconstruction artifacts.

Discussion and conclusion
In this work we provide the proof of concept that a combined parametric and non-parametric approach is feasible for accurate reconstructions from a single in-line hologram. The parametric part of the reconstruction fits the high contrast fringes thereby enabling the non-parametric part to retrieve the faintest parts of the object. These reconstructions are performed jointly, meaning that their interferences are taken into account, in contrast to the previous method that works on subtracted intensities [9].
The choices of the Mie model and the Rayleigh-Sommerfeld propagation make the proposed method adapted to a large domain of validity extending to absorbing and/or dephasing spherical objects and transmittance planes without any restriction on the recording distance. Thus, other applications can be considered such as lens-free microscopy [4,5,7,34] or in-line microscopy in biology [11,14,24] where the Lorenz-Mie model fitting and the non-parametric reconstruction techniques have independently shown their efficiency. The proposed method would be particularly well suited to reconstruct samples in which spherical objects are present among objects of more complex shapes [35,36].
Finally, some work remains, especially on the non-parametric part of the algorithm. Indeed, some signal remains in the residues, in high intensity area such as close to the droplets or in the highest concentration parts of the vapor wakes. This suggests that a better choice of regularization and associated hyper-parameters would improve the reconstruction to more accurate quantitative values.

Appendix A: Forward model
This appendix details the forward model used in section 3. Placed in the context of monochromatic plane wave illumination of wavelength λ, it can be assumed without loss of generality that it is normalized to U inc = 1 on the sensor plane z = 0. Thus, at a given geometrical distance z > 0 of this plane, the incident wavefront is: where k = 2π/λ is the wavenumber and the minus sign accounts for the propagation of the illumination toward the sensor.

Parametric step: model fitting of the Lorenz-Mie theory
In our application in fluid mechanics, the parametric part u of the object is a set of droplets of ether. Assuming they are homogeneous spheres with a known complex refractive index, their diffraction can be modeled via the Lorenz-Mie theory [19,[37][38][39][40]. This provides an analytic solution for the diffracted electromagnetic fields E di f and H di f scattered by a sphere, in spherical coordinates (r, θ, ϕ) centered on the sphere center. In the case of plane wave illumination, each component of the electromagnetic field is written as finite series over a single integer n, involving Ricatti-Bessel function ξ n , Legendre polynomials π n and τ n and the Mie scattering coefficients a n and b n which depend only on the wavelength λ, the sphere radius and its refractive index. For example, the θ component of the scattered electric field E θ di f (r, θ, ϕ) is written as: pw n a n ξ n (kr) τ n (cos θ) cos ϕ − ib n ξ n (kr) π n (cos θ) sin ϕ , where E 0 is a normalization factor and c pw n = (−i) n 2n+1 n(n+1) . A change of variable provides the electric and magnetic field components in a Cartesian coordinates system (x, y, z).
In the present case, the distances between the droplets and the sensor are large. As a consequence, the problem can be expressed in the realm of the scalar theory of light propagation [15]. In that case, the complex scalar amplitude can be identified with a given component of the Cartesian vector field, for instance the one on the x-axis: U Mie di f = E x di f . Thus, the light diffracted by N d spherical droplets is described by a limited number of parameters: respectively their 2D position on the sensor plane, their distance to the sensor, their radius and their complex refractive index. In the present situation, the diffraction pattern does not depend on the refractive index due to the long distance between the droplets and the sensor. The light diffracted by the droplets is lost and as a consequence they behave as opaque objects. Thus their refractive index n d is removed from the set of unknown parameters.
For each particle in the field, the parameters u d are locally optimized with the Lorenz-Mie model, in an iterative inverse problem approach with a numerical estimation of the gradient at each step. Note that even if a simpler formalism would be sufficient in the present case [13,20], the Lorenz-Mie model is kept, without far field approximation, to provide the most general framework. Thus: and at a given iteration of the alternative approach 1, u is updated as follows: As mentioned above, this kind of parametric model fitting is robust. Thus, this optimization (10) is not sensitive to the minimization on c in the data fidelity term definition (3). In addition, in the present case, most of the signal in the data consists of high contrast fringes produced by the droplets. These fringes being efficiently "explained" by the parametric step, the data is correctly fitted in terms of contrast. Therefore the estimation of the parameter c in the non-parametric step is better constrained.
The proposed method consequently provides a way to calibrate the acquisition dynamics.
Non-parametric step: Rayleigh-Sommerfeld propagation for an image-based reconstruction

Diffraction model
In the present application in fluid mechanics, the non-parametric part v of the object is linked to the vapor wake generated by evaporation of the ether droplets. Then, an intuitive unknown would be a 2D map of its complex transmittance t 2D . Under the hypothesis of an optically thin medium, it represents the projection of the sample along the line of sight. The diffracted wave U np di f is produced by the difference in the transmittance to the unitary plane: t 2D − 1. After propagation from the object plane at a distance z, the resulting diffracted wave on the sensor plane is given by the convolution: where the convolution kernel h z (x, y) is the Rayleigh-Sommerfeld propagator [15] at the illumination wavelength λ and for a propagation distance z: where r = x 2 + y 2 + z 2 .
To describe the complex transmittance plane, by splitting the absorbing and dephasing properties of the object, t 2D is actually modeled by its complex optical length difference : t 2D = e 2i π λ . The real part of consequently defines the dephasing properties of the object, and its imaginary part the absorption.

Constraints and regularizations
As mentioned earlier, minimizing (4) implies adding adequate prior information about the object.
In the present case, the droplets are the only objects that behave as absorbing ones, described by u in the parametric part of model. Thus, their vapor wakes and the surrounding medium can be considered as purely phase objects, described by a real-valued difference in optical length, leading to the following constraint on the imaginary part of : In terms of regularization, the assumption of an optically thin medium supposes a slow spatially varying object v. This is favored by a low energy gradient. For a given 2D real-valued image M, it consists in minimizing the following regularization term: where the directional gradients on x ↔ j and y ↔ i are computed by finite differences between adjacent pixels and where W is a weighting matrix. This assumption is nevertheless not fulfilled in the vicinity of the absorbing droplets that introduce a high contrast change in the transmittance plane co-localized with their edges as can be seen in Figs. 2(b)-2(d) and Figs. 3(c)-3(n). In other words, the opaque droplets mask the vapor wake, which induces a zero-valued phase on their support, i.e. a sharp edge at the frontier, as shown in Fig. 2(d). Favoring sharp-edged objects is achieved by concentrating the energy of its gradient. For a given 2D real image M, it consists in minimizing the following regularization term, to enforce a sparse gradient [41,42]: where the directional gradients on x ↔ j and y ↔ i are computed by finite differences between adjacent pixels, W is a weighting matrix and is a small real-valued constant to relax the regularization and insure its differentiability. An adequate trade-off is required to combine these regularizations. To this end, the parameters u estimated by the step (10) are used to build specific weighting matrices W d r for the N d droplets acting as a mask on the droplets, as presented in Figs. 5(e)-5(f): • W d r = 1 on a disk of radius η in r d around the d th droplet of radius r d , • W d r = 0 outside a disk of radius η out r d around the d th droplet of radius r d , • a power law is used to obtain a smooth transition between the inner disk and the outer disk.

Multi-z reconstructions
A final ingredient can be added to the model to correctly account for the strong diversity of the droplet distances. Due to the low spatial frequencies of the vapor wakes in the context of the optically thin medium, the propagation distance z is not critical for the range of droplet depths. Nonetheless, the high spatial frequencies, i.e. sharp edges, held in the vicinity of the droplet make this reconstruction distance critical to correctly focus the shadowod of the droplets in the transmission plane. Figure 4 highlights the need to correctly account for the different propagation distances. In the case where the reconstruction plane is set at an unique droplet distance, like in Figs. 4(b),4(d),4(g),4(h), only the corresponding droplet shadow is correctly reconstructed. The high frequencies produced by the other droplets cannot be explained by the model: the residues worsen and a defocused image of the droplets shadows impacts the reconstruction. Conversely, with the method proposed in the following, the different propagation distances are taken into consideration. In so doing, Figs. 4(a),4(c),4(e),4(f) show that the phase map is discontinuous as expected at the droplet positions and the high contrast fringes are better explained in the residues.
Thus, if several droplets are present at different distances in the field, the forward model (11) must be adapted accordingly. To do so, N d patches are built for the N d droplets, having the same support as their corresponding weight W d r . Under the first Born approximation [30], each of these patches d scatters the unperturbed incident wavefront U inc (z), assuming that the latter was not modified by the upstream planes. All the resulting waves interfere on the sensor plane, as sketched in Figs. 5(a)-5(d). By convention, the slowly evolving phase is placed on the plane z 1 of the first droplet. According to the optically thin medium hypothesis, the overall optical length is = N d d=1 d .
2π λ d . (e) Weighting matrix N d d=1 W d r that makes it possible to combine edgepreserving regularization around the droplets (w = 1, white) and smoothing regularization on the background (w = 0, black). In the transition areas, the two regularizations are combined. (f) Profile of the weighting coefficients drawn along the dashed line of (e) on the first droplet. The droplet radius r 1 is emphasized as are the boundaries defined by η in r 1 and η out r 1 , linked by a power law.
Thus, defining the non-parametric part of the object v as the set of all the phase shifts introduced by the optical lengths d : v = {v d = 2π λ d } d ∈ 1, N d , the diffracted wave is: where z d are the distances obtained in the parametric step.
Combining the constraint (13), the regularizations (14) and (15) and the propagation model (16), at a given iteration of algorithm 1, v is updated as follows: with the constraints D on the domain of v: and where µ 1 and µ 2 are hyperparameters to balance the regularization terms compared to the data fidelity term. It should be mentioned here that for any v, an analytic solution for the scaling factor c in the definition of the data fidelity term (3) can be computed on the fly. Indeed, the optimal scaling factorĉ is given by: The minimization problem (17) is solved with a limited-memory quasi-Newton method with bound constraints, VMLM-B [43]. We use the version implemented in the GlobalBioIm framework [44]. Solving (17) implies evaluating the gradient of the data fidelity term. This term contains a squared complex modulus depending on the complex-valued unknowns. The derivation of the close-form expression of the gradient is detailed in appendix.

Numerical reconstructions
In this paper, all the reconstructions presented are run with the following empirically chosen parameters: • a first rough estimation u 0 of the parametric part u is performed following a matching pursuit algorithm by browsing the discretized space of the parameters around the radii r d and the distances z d estimated via a simple backpropagation of the hologram, • v is initialized to null dephasing planes v 0 = {v d = 0} d ∈ 1, N d • c is estimated on the fly via Eq. (19), • η in = 1.5 and η out = 12.5 with a second degree power law to link the two regularizations, • At each loop, N np it = 50 iterations are performed to solve the problem (17), • In equation 11, the propagation -a discrete convolution by the kernel h z -is performed in the Fourier space with fast Fourier transforms involving a zero-padding operation to avoid aliasing effects. This final operation can also be used as a masking step in the forward model to allow reconstruction of the vapor wakes on a greater field of view than the one seen by the sensor [13,31].