Accurate 3D tracking and size measurement of evaporating droplets using in-line digital holography and “inverse problems” reconstruction approach

Mozhdeh Seifi; Corinne Fournier; Nathalie Grosjean; Loic Méès; Jean-Louis Marié; Loic Denis

doi:10.1364/OE.21.027964

1. Introduction

In fluid mechanics, the study of bubbles and droplets carried by flows essentially relies on velocity and size measurements. The most widely used measurement technique in this context is Phase Doppler Anemometry (PDA) [1]. PDA provides only single-point measurements and therefore is incompatible with the tracking of individual particles (i.e., Lagrangian tracking). Interferometric Laser Imaging for Droplet Sizing (ILIDS) [2] and Global Phase Doppler Anemometry [3] provide alternative solutions to measure droplet (or bubble) size and location in a whole flow section. In principle, these techniques can be extended to estimate 3D location of particles, but in practice their depth of field is limited by the small thickness of a laser sheet. Small displacements in depth (i.e., parallel to the optical axis) can be considered in order to measure a third component of velocity [4]. Such measurements still remain essentially 2D as the measurement volume is flat, with a field of view much wider along the transversal directions than in depth. Lagrangian tracking of evaporating droplets requires time-resolved 3D imaging system with a full 3D measurement volume, large enough to contain a significant part of a droplet trajectory. To this purpose, digital holography is a very promising technique, allowing both 3D location and droplet size tracking with a good temporal resolution and a wide/deep field of view.

Digital holography is being increasingly used in applications that require micro-objects tracking (e.g., [5, 6, 7] ). In contrast to off-axis setups, the in-line setup (i.e., the Gabor setup) is less sensitive to vibrations because it does not involve beam splitters, mirrors or lenses. It also exploits the whole frequency bandwidth of the sensor to encode accurately the depth and the size of the objects on the holograms. This imaging technique is also called “lensless imaging” [8, 9, 10] , as it involves no lens between the object and the sensor. A key issue is to design hologram processing algorithms that achieve the best possible accuracy in 3D location and sizing. Over the past decade, numerous algorithms for the analysis of digital holograms have been proposed and several journal special issues were published on the subject, e.g., [11, 12, 13]. The most widely used approach is based on the simulation of the optical reconstruction of the hologram (using, e.g., Fresnel functions to back-propagate the hologram), following by segmentation of this 3D image and estimation of the object parameters. These techniques have a limited accuracy due to the following characteristics of the imaging technique [14, 15] : (i) the signal is truncated resulting in biased estimations, especially for objects located at the borders of the field, (ii) the spatial resolution of the sensor is low which either results in the presence of ghost images in the reconstruction or forces the experimenter to increase the working distance which reduces the hologram signal to noise ratio and thus degrades the accuracy of estimation. Furthermore, most of these reconstruction algorithms require the user to tune several parameters. These methods already achieve very interesting results [16, 17, 9, 18].

Another family of hologram processing methods follow a signal processing approach to perform the detection and sizing of the objects directly from the hologram rather than from the back-propagated optical field [19, 20, 21, 22]. In this article, an image processing algorithm based on the “inverse problems” approach [23, 24] is used to perform an accurate estimation of both size and 3D position of evaporating diethyl ether droplets. As the Gabor holographic setup is restricted to small objects at low concentration (satisfying Royer criterion [25]), the hologram can be approximated as the sum of the diffraction patterns produced by the objects. The “inverse problems” methods aim to invert the hologram formation finding a set of diffraction patterns that closely fit the experimental data. In the case of spherical objects, finding the best fit corresponds to estimating the size and 3D position parameters. In contrast to the light back-propagation approaches, these calculated patterns intrinsically take into account signal truncation and the low spatial resolution of the sensor. In addition to an improved accuracy of estimation, such unsupervised approaches can greatly expand the field of view outside of the sensor area [24, 26].

In this paper we employ a previously introduced parameter estimation method [23, 24, 21] , to study the evaporation of diethyl ether droplets from hologram videos. The droplets are assumed to be spherical. The parametric description of each object then relies on four parameters: three spatial coordinates and a single shape parameter (the radius). We chose ether to evaluate the performance of the method, because it evaporates very fast.

The structure of this paper is as follows: the recording setup is detailed in Sec. 2. In Sec. 3 the image formation model and the iterative model fitting algorithm are presented as well as two improvements of the algorithm regarding accuracy and time costs. The reconstruction of experimental holograms, presented in Sec. 4, shows that in-line digital holography and the reconstruction algorithm are suitable for the study of evaporation processes. The paper is concluded in Sec. 5.

2. Description of the experimental setup

An in-line digital holography setup, also called the Gabor setup, is used in this paper to record the holograms. The diethyl ether droplets are generated by a piezoelectric jetting device manufactured by MicroFab Technologies (see Fig. 1). This injector generates mono-dispersed droplets with radii of 31 μm ± 0.25 μm [21]. The droplets are injected at distances of 45 cm to 52 cm from a 800 × 1280 pixel Phantom V611 camera with the pixel size of 20 μm and the fill-factor (i.e., active area over total area of pixel) of 0.56. The frame rate is set to 620 frames per second. The illuminating laser beam is produced using a Nd:YVO4 laser (Spectra-Physics, Millenia). The laser beam divergence introduces a magnification factor of approximately 2.4 in the system (see Sec. 4.1). The experimental holograms have a signal to noise ratio (i.e., the ratio of the magnitude of the signal to the standard deviation of noise) ranging from 5 to 9. Given the magnification factor and in order to investigate the whole evaporation process, it proved necessary to capture three sets of holograms to investigate the whole evaporation process. The set “0” captures the holograms of the droplets being injected into the air. The injector’s diffraction pattern is visible on this set. The sets “1” and “2” capture the holograms with relative camera translations of 7.5 mm and 15 mm compared to the first setup (see Fig. 2). For each set, a video of 100 holograms is recorded. One hologram of each set after removing its background is shown in Fig. 2.

Fig. 1 A picture of the droplet jet

Download Full Size | PDF

Fig. 2 The setup for capturing three sets of holograms containing evaporating diethyl ether droplets. The droplets radii are scaled up for the sake of visualization.

Download Full Size | PDF

3. Automated estimation of particles parameters

Signal processing methods can be employed in the case of parametric objects like droplets to estimate the size and 3D position of every droplet directly from the hologram. In the case of low density object fields (i.e., satisfying the Royer criterion [25] ), the image formation model can be considered as linear and the hologram can be decomposed into the sum of the diffraction patterns of the objects [23]. When objects are modeled with a parametric shape, the model of their diffraction pattern is parametric, too. It depends on the shape and 3D position of the object as well as on the experimental setup parameters. A previously introduced iterative algorithm [23, 24] finds the size and the coordinates of the droplets with high accuracy.

In the next subsections, the hologram formation model of a small spherical opaque object is presented and a background removal method is deduced from this model. The iterative algorithm is then briefly revisited and two improvements are detailed which make the parameter estimation more accurate and faster.

3.1. The parametric diffraction pattern model

We begin by describing the forward model of hologram formation in order to build the parametric model of diffraction patterns. The diffraction pattern of a spherical droplet of radius r_n located at distance z _n from the sensor is here assumed to be identical to the diffraction pattern of a circular opaque disk of the same radius at the same distance from the sensor. Let us define the aperture of the opaque disk ϑ_n as equal to 1 inside a disk (with radius r_n) centered on 2D position (x_e, y_e) and 0 outside this disk. Under the far field condition (4 πr_n²/(λz_n) ≪ 1), the complex amplitude of the signal A_holo (x, y) on every pixel (x, y) of the camera can be written as [27]:

A_{holo} (x, y) \approx A_{ref}^{0} (x, y) - A_{ref}^{z_{n}} (x_{n}, y_{n}) \cdot (ϑ_{n} * h_{z_{n}}) (x, y),

where

A_{ref}^{0} (x, y)

stands for the complex amplitude of the reference wave on the hologram plane and

A_{ref}^{z_{n}} (x_{n}, y_{n})

refers to the complex amplitude of the reference wave on the object aperture plane and is assumed to be uniform on the object’s aperture. In this formulation, *corresponds to the two dimensional convolution operator. Besides, h _{z_n} (x, y) represents the Fresnel function:

h_{z_{n}} (x, y) = [1 / (i λ z_{n})] exp [i π (x^{2} + y^{2}) / (λ z_{n})],

where λ stands for the laser wavelength and

i = \sqrt{- 1}

.

The captured intensity on the sensor I _holo (x, y) is calculated as

I_{holo} (x, y) = {| A_{holo} (x, y) |}^{2},

leading to:

I_{holo} (x, y) = I_{ref}^{0} (x, y) - 2 ℜ {A_{ref}^{0 *} (x, y) \cdot A_{ref}^{z_{n}} (x_{n}, y_{n}) \cdot [ϑ_{n} * h_{z_{n}}] (x, y)} + β .

In this expression β refers to the second order term which can be considered negligible for small particles with $r_{n} ≪ \sqrt{λ z_{n}}$ .

In the case of a spherical diverging reference beam with a point source located on the optical axis at a distance z _s from the sensor, this expression can be simplified (see Appendix. A ):

I_{holo} (x, y) \approx I_{ref}^{0} (x, y) - \sqrt{I_{ref}^{0} (x, y)} α_{e} . [ϑ_{e} * ℜ (h_{z_{e}})] (x, y) .

In this equation α_e is a proportionality factor accounting for the non-uniformity of the illuminating beam and z_e = mz_n, where $m = \frac{z_{s}}{z_{s} - z_{n}}$ is called the magnification of the system and can be estimated through a calibration procedure (see Sec. 4.1).

A spherical droplet with an aperture ϑ_n is therefore reconstructed as an object of aperture ϑ_e with a radius of r_e = mr_n located at x_e = mx _n , y _e = my _n and z _e. Defining the set of parameters θ _e = ( x _e , y _e , z _e , r _e ) and the oscillating term g _{θ_e} = −[ ϑ _e ℜ (h _{z _e} )]( x , y ), the intensity expression is simplified to:

I_{holo} (x, y) = I_{ref}^{0} (x, y) + α_{e} \sqrt{I_{ref}^{0} (x, y)} \cdot g_{θ_{e}} (x, y) .

Under the assumption of having $r_{e} ≪ \sqrt{λ z_{e}}$ , the diffraction model can be mathematically formulated as [28] :

g_{θ_{e}} (x, y) = - \frac{π r_{e}^{2}}{λ z_{e}} sin (\frac{π ρ_{e}^{2}}{λ z_{e}}) J_{1 c} (\frac{2 π r_{e} ρ_{e}}{λ z_{e}}) {Pix}_{θ_{e}} (x, y),

with

ρ_{e} = \sqrt{{(x - x_{e})}^{2} + {(y - y_{e})}^{2}}

. In this formulation, Pix _{θ_e} ( x , y ) accounts for the pixel integration performed by the sensor on the active area of the pixel:

{Pix}_{θ_{e}} (x, y) = sinc (\frac{κ π (x - x_{e})}{λ z_{e}}) sinc (\frac{κ π (y - y_{e})}{λ z_{e}}),

where κ denotes the width of the square active area.

3.2. Non-uniform background removal

In practice, having holograms with non-uniform backgrounds is common. Such a background image can degrade the accuracy of parameter estimation by introducing correlated noise. Nonuniform background can be produced by non-uniform laser beams, scratches or dust on the cover glass of the camera or environmental noise. A non-uniform background is classically removed either by subtraction or division. This background image is either captured in the absence of object or is estimated as the mean image of a video of holograms.

An effective background removal method can be derived from the image formation model of Eq. (5) : the background image I _ref is first subtracted from the hologram and the resulting intensity is divided by the amplitude of the background, leading to the pre-processed hologram d ( x , y ):

d (x, y) = \frac{I_{holo} (x, y) - I_{ref}^{0} (x, y)}{\sqrt{I_{ref}^{0} (x, y)}} = α_{e} . g_{θ} (x, y) .

For each measurement set in this paper, a distinct background image is calculated corresponding to the mean of 100 holograms. An example of a droplet hologram after background removal is provided in Fig. 3(a). Such an operation leads to holograms with uniform backgrounds if objects are distributed randomly in the volume. In our case, a dark mark can be noticed along the trajectory of the droplets, because droplets locations do no cover uniformly the sensor during the time sequence, but rather are aligned along the same trajectory (see Fig. 3(a-1) ).

Fig. 3 Illustration of a pixel region exclusion mask application on a hologram of set “1”. (a-1) an experimental hologram, (a-2) the experimental hologram being masked, (b-1) cleaned hologram, (b-2) masked cleaned hologram. The colored pixels show the masked pixels which were not included in the parameter estimation.

Download Full Size | PDF

3.3. Iterative parameter estimation algorithm

In the case of volumes with a low density of objects, an iterative algorithm [23, 24] can be used to estimate the parameters of the diffraction pattern of each object. This algorithm performs the parameter estimation task, object by object, by minimizing the weighted least square difference between the experimental data and a parametric model. The method is optimal from the signal processing point of view if noise can be considered as white and Gaussian. Having such assumptions on noise, we introduce a binary mask w to exclude, from the analysis, some regions of the sensor that do not contain meaningful data (see Sec. 3.4 for a detailed discussion on the mask). The binary mask can be generalized to a weighting mask to account for non-uniformity of the noise variance in the image, as described in [24 ].

To simplify the notation, we define the masked zero-mean image (denoted by overlining the image variable u ) ū

\bar{u} (x, y) = u (x, y) - \frac{1}{c} \sum_{x} \sum_{y} w (x, y) u (x, y),

where c represents the number of pixels included in the mask:

c = \sum_{x} \sum_{y} w (x, y) .

Using this notation, the mentioned iterative algorithm consists of three main steps:

an exhaustive parameter estimation step which searches in a sampled parameter space { θ = ( x , y , z , r )} to find a rough estimation of the parameters of the current particle. It has been previously shown [23, 24] that minimizing the least squares difference is equivalent to maximizing the weighted normalized correlation term between the data and the model. The parameter set θ ̂ _e that results in the maximum correlation (denoted by arg max in the following) is found from the following constrained optimization problem: ${\hat{θ}}_{e} = \underset{θ = (x, y, z, r)}{arg max} {\frac{{[\sum_{x} \sum_{y} w (x, y) \bar{d} (x, y) {\bar{g}}_{θ} (x, y)]}^{2}}{\sum_{x} \sum_{y} w (x, y) {\bar{g}}_{θ}^{2} (x, y)}}, s . t . α \geq 0,$ where $α = c \frac{\sum_{x} \sum_{y} w (x, y) \bar{d} (x, y) {\bar{g}}_{θ} (x, y)}{\sum_{x} \sum_{y} w (x, y) {\bar{g}}_{θ}^{2} (x, y)} .$ In this formulation, θ̂ _e corresponds to the rough estimation of the object parameter θ _e. A diffraction pattern with inverted contrast should not be detected as an object signature. Such patterns are excluded from consideration by enforcing α ≥ 0.
a local optimization step which aims to refine the estimated parameters of the previously found pattern by slightly changing the parameters in a continuous parameter space. This step finds the parameter set that minimizes the following nonlinear weighted least squares problem (denoted with arg min) in the neighborhood of θ̂ _e : $θ_{e}^{*} = \underset{θ = (x, y, z, r)}{arg min} {\frac{1}{c} \sum_{x} \sum_{y} w (x, y) {[\bar{d} (x, y) - α_{e}^{*} {\bar{g}}_{θ} (x, y)]}^{2}}, s . t . α_{e}^{*} \geq 0.$ In this formulation, $θ_{e}^{*}$ corresponds to the final estimate of θ _e. The rough estimation of the parameters of the first step θ̂ _e must be accurate enough to be in the main basin of this non-convex cost function. During the optimization process, $α_{e}^{*}$ is updated at each iteration according to Eq. (13).
a cleaning step which consists in subtracting the weighted signature of the found object (i.e., $α_{e}^{*} {\bar{g}}_{θ_{e}^{*}}$ ) from the hologram in order to increase the signal to noise ratio of the remaining objects. The algorithm then starts over to find the next object from the residuals.

Let us note that a last refinement step can be performed by cleaning all the detected particles except one, and then carrying out step (ii) in order to refine the parameters of this particle with an improved signal to noise ratio.

To speed up the weighted normalized correlation computations (of step (i)), fast Fourier transforms are used [23, 24]. A free Matlab ^® toolbox called “HoloRec3D” implements the iterative algorithm for digital hologram reconstruction of spherical micro-particles [29, 30].

3.4. Pixel region exclusion mask

The first advantage of the binary mask w used in step (i) (see Eq. (12) and Eq. (14) ) is to take into account the truncation of the data using FFTs for faster calculations. The second advantage is to exclude pixel regions [in steps (i) and (ii)] where the parametric model can’t explain the data. A closer look to Eq. (7) reveals that the intensity of the model is equal to zero at the center of the pattern (i.e., x = x _e and y = y _e ). However the experimental holograms of evaporating droplets show high intensities in the pattern centers (see Fig. 3(a-1) ) which may be caused by a film of diethyl ether vapor surrounding the droplet during evaporation. The diffraction pattern produced by the vapor film is not included in the diffraction model used here. Because it is spatially well localized, this effect can be removed from consideration by using a binary weighting mask. The binary mask contains the values ”one” for the useful pixels, and ”zero” for the pixels in a circular neighborhood of the center of evaporating patterns, in order to discard regions where our model is too coarse (see Fig. 3(a-2) for a masked hologram). The fit between the experimental data and the model is thus performed only on the high frequencies of the pattern.

To create the mask, a first rough estimation of the parameters is performed to obtain the centers of the masking areas (i.e., ( x _n , y _n )). The size of the circular neighborhood is the same for all diffraction patterns. The best size is the one that gives the highest weighted normalized correlation value between a test hologram and the fitted model. One residual hologram from which the estimated models have been subtracted (during the cleaning step), is shown in Fig. 3(b-1) . The high frequency part of the diffraction patterns is well cleaned. It illustrates that the model accurately fits the data (see also Fig. 3(b-2) ) in high frequency regions. Radial mean profiles of experimental data and the diffraction pattern model are presented in Fig. 4 together with the residuals.

Fig. 4 Comparison of the radial mean profile of one evaporating droplet’s signature (calculated from a hologram of set “2”) with the fitted model, the residuals and the masked part of the signal. (b) shows a zoomed-in version of (a).

Download Full Size | PDF

The exclusion mask allows parameter estimation for evaporating droplets, using a simple hologram model. As a future work, another approach would consist in using a more appropriate model, taking into account the effect of the surrounding vapor film. Such a model has been very recently proposed by our group [31]. Using this new model would require the estimation of more parameters, which would be computationally more demanding, yet it could lead to a richer description of the evaporation process.

3.5. Rough parameter prediction

The most costly step of the iterative algorithm described in Sec. 3.3 is the exhaustive parameter estimation step (i) as it searches in the whole parameter space to find a first and rough estimation of the parameters (e.g., for ( x , y ) parameters, the search is carried at least on the whole sensor). Such a roughly estimated parameter set is next used as an initial guess for the local optimization step (ii). Recently, two algorithms have been proposed to reduce the time costs of the exhaustive search [28, 32] based on multi-scale approach or using low-dimensional approximation of diffraction pattern models. In the case of videos of droplet holograms, a rough physical model of the trajectories can be used to coarsely anticipate the parameters of the objects at time t from the previously estimated parameters. An unsupervised heuristic anticipation routine is introduced in this section. Considering the object parameters as θ , the goal is to predict the size and location of all objects (denoted as θ̃ _t for each object) using the estimates computed from the previous holograms i.e., $θ_{t - Δ t}^{*}$ (Δ t stands for the acquisition period). As diethyl ether droplets are in free fall, a rough estimation of the parameters of a droplet at time t can be predicted by using a second degree polynomial model describing the motion of the droplets. The prediction of the whole parameter vector is thus modeled using:

{\tilde{θ}}_{t} = θ_{t - Δ t}^{*} + \frac{1}{2} a_{t - Δ t} Δ t^{2} + v_{t - Δ t} Δ t,

in which the ”velocity” and ”acceleration” vectors, v _t and a _t , have four components (three position coordinates and a radius), deduced from the previous accurate coordinates

θ_{t}^{*}

,

θ_{t - Δ t}^{*}

and

θ_{t - 2 Δ t}^{*}

that are estimated by the iterative algorithm:

a_{t} = \frac{θ_{t}^{*} + θ_{t - 2 Δ t}^{*} - 2 θ_{t - Δ t}^{*}}{Δ t^{2}},

v_{t} = \frac{θ_{t}^{*} - θ_{t - Δ t}^{*}}{Δ t} .

The first two holograms of each set are processed using the method described in Sec. 3.3 to provide the initial information for the prediction of the parameters for the next holograms. When the anticipation is not accurate enough to get an initial point in the main basin of the cost function (i.e., the optimization does not converge), this anticipation step is discarded and the exhaustive search (i) takes over to find the initial point for the optimization. Our results however showed that for 90% of the cases, the anticipation formula was accurate enough to perform a successful optimization step. This simple prediction method reduced the time costs of accurate parameter estimation from 6 minutes to 30 seconds.

4. Experimental results

In this section, the results of volume reconstruction using the iterative algorithm are presented. First, in Sec. 4.1, the calibration process is described. In Sec. 4.2 water droplet trajectories are reconstructed and size measurements are compared to PDA measurements. In Sec. 4.3 the iterative algorithm is applied to evaporating ether droplets, showing that the technique provides results that are accurate enough to study the evaporation phenomenon (e.g., the evaporation rate).

4.1. Calibration

To perform the calibration step, holograms of a glass reticle with a linear scale (Edmunds Optics, #62–252) are captured with varying distance between the reticle and the sensor. The reticle is first placed at a distance of z ₀ = 473.5 mm from the sensor and moved back and forth in a depth range of 10 mm by steps of 1 mm.

For each position of the reticle, the hologram is re-focused numerically using a standard reconstruction algorithm. An example of the in-focus reticle image is shown in Fig. 5(a). The reticle’s grading is next used to calculate the magnification with an estimated uncertainty of 0.005. A regression routine is finally used to calculate the magnification as a function of depth as shown in Fig. 5(b). The regression results in the following magnification formula which is used later to transform the estimated parameters ( x _e , y _e , z _e , r _e ) into the object reference frame coordinates ( x _n , y _n , z _n , r _n ):

m = 7 (z_{n} - 473.5 \times 10^{- 3}) + 2.36 .

Fig. 5 Calibration process, (a) a sample of a reconstructed calibration hologram, (b) the fitted regression curve used to estimate the magnification.

Download Full Size | PDF

4.2. Non-evaporating droplets

PDA can be considered as a reference method for spherical particle size measurement. We use it to validate the accuracy of size estimations obtained by our digital holography method. PDA measurements are performed at the outlet of the injector just before recording the holograms. The mean droplet diameter, estimated by PDA from 3500 measurements, is 31.07 μm ± 0.078 μm. These values are close to those found in the same conditions in another work [ 21] , showing that the injection is well reproducible. Using the setup described in Sec. 2, 310 droplets are reconstructed. The mean radius and the standard deviation of the radius are estimated as 30.57 μm and 0.11 μm, respectively. This proves that the accuracy of digital holography is comparable to that of PDA, with a relative standard deviation of 0.3%. One of the holograms of such droplets is shown in Fig. 6(a). Figure 6(b) shows the residuals after cleaning the diffraction pattern of the droplet.

Fig. 6 A hologram of non-evaporating droplets (a) and the residuals after cleaning the diffraction pattern of the droplet (b).

Download Full Size | PDF

4.3. Evaporating droplets

The three sets of holograms described in Sec. 2 are processed using the iterative algorithm (Sec. 3.3). The droplets of set “0” are just leaving the outlet of the injector and their velocity is high, typically 0.23 m.s ⁻¹. As a consequence, the reconstructed trajectories from set “0” consist of only 5 or 6 droplets positions and are therefore not discussed here. As previously mentioned in Sec. 3.5, the first two holograms of every set were processed using the exhaustive search step. The x,y dimensions of the search space were sampled by steps equal to the pixel size. In our experiments, z and r dimensions were sampled by steps of 0.8 mm and 1.2 μm respectively. Estimated parameters of the first two holograms of each set were used to compute rough estimations of the parameters of the next holograms, using the anticipation method. The algorithm therefore skips the costly exhaustive search step (see Sec. 3.5).

The trajectories containing more than 10 droplets are reconstructed from set “1” (17 trajectories with an average of 34 droplet positions) and set “2” (13 trajectories with an average of 16 droplet positions). 3D visualization of the trajectories for set “1” and set “2” are shown using videos in Figs. 7(a)–7(b) respectively. Some examples of the squared radius evolution over time are presented in Figs. 8(a)–8(b) for set “1” and set “2” respectively. We note that after some time, the radius stabilizes around a constant value, on the order of 10 micrometers for all the runs (see Figs. 8(a)–8(b) ). We suppose that the evaporation stops at this stage, when all the ether has evaporated. The remaining droplet would then be composed of only water and would not evaporate anymore. The percentage of water given by the manufacturer, and experimentally verified, is 0.2%, while based on the final radius estimation, an initial concentration of about 4% is measured. Because these experiments are performed in humid air (with relative humidity of 31.8%), a plausible explanation for the higher percentage of water is that the fast evaporation of the diethyl ether cools the humid air around the droplet, thus causing some condensation at the surface. The existence of such condensation phenomenon was reported by Law, [33 ] for alcohol droplets vaporizing in humid air. The repeatability of the droplet size evolution over time for various 3D trajectories can be noted from Figs. 8(a)–8(b).

Fig. 7 3D visualization of reconstructed trajectories (a) movie for set “1” ( Media 1 ), (b) movie for set “2” ( Media 2 ). Scale of the sphere radius has been increased to visualize the evaporation phenomenon.

Download Full Size | PDF

Fig. 8 Evolution of the radius squared of the droplets over time (a) set “1”, (b) set “2”

Download Full Size | PDF

The droplets that contain only water (i.e., the ones with constant radius) can be seen in Fig. 3(a-1). They correspond to the four patterns located at the bottom, which do not exhibit the high contrasted central disturbances created by the vapor film, confirming our belief that evaporation has stopped. The use of the mask is not necessary for such water droplets and we note that they are well cleaned from the residual image (see Fig. 3(b-1) ). Further investigation could be done to confirm this assumption, by measuring the refractive index of the droplets along their trajectories or by varying the humidity and temperature of ambient room conditions.

The droplet trajectories are then analyzed and results are compared to a simple evaporation model. When evaporation has reached stationary conditions (i.e., reaching constant droplet temperature), the squared diameter of a droplet generally decreases linearly with time. In terms of squared radius, this so-called “ d ² ” law [34] is expressed by

r {(t)}^{2} = r {(0)}^{2} + K t,

where K < 0 is the evaporation rate. This linear decrease can be observed for big evaporating droplets in Fig. 8 . A theoretical value of K can be calculated as a function of the physical properties of diethyl ether, temperature and pressure. As the temperature of the droplets at the output of the injector is difficult to estimate, we can only give a range for K values: [−7.5, −6.2] μm ² / ms (see Appendix B for the detailed calculation of this range).

To estimate K from the results obtained from digital holography, a least squares linear fit was performed on the first part of the squared-radius-over-time curves (i.e., the part of the curves that presents evaporation). Figure 8 shows the regression lines which closely fit the data. The average K is found equal to −7.7 μm ² / ms and −6.2 μm ² / ms for set “1” and set “2” respectively.

To calculate the accuracy of the measurements, each trajectory is used to compute the distances between the estimated radii and the fitted curve [35 ]. The standard deviation of these distances gives a rough estimation of the accuracy on the radius measurement. For set “1”, the standard deviation is 0.5 μm for evaporating droplets with the radii in the range [9, 26] μm. The same study on set “2” gives a standard deviation of 0.4 μm in a radius range [13, 18 ] μm.

This evaluation of the precision overestimates the actual precision, since all deviations from the evaporation model are attributed to measurement errors [35]. This precision is still much larger than the precision estimated with non-evaporating droplets. The most reasonable explanation is that our model of the diffraction pattern of an evaporating droplet is coarse. The exclusion mask w accounts for the large misfit of the model at the center of the diffraction pattern. The higher frequency part of the model (i.e., corresponding to the larger diffraction rings) is also distorted, probably because of the presence of vapor, leading to a poor fit in this region [see Fig. 3(b-2)], and to coarser estimates of the depth and radius parameters. By excluding some pixels from the estimation procedure (about 7% of the pixels), the mask also degrades the accuracy.

5. Conclusion

Fast evaporation of free-falling diethyl ether droplets has been studied by means of in-line digital holography. To perform the reconstruction, an iterative model fitting algorithm based on an “inverse problems” approach [23, 24] has been used. This algorithm has been already proven to result in accurate 3D reconstructions of particles. This approach requires a mathematical parametric model of objects’ holographic patterns. In the case of evaporating droplets, the presence of vapor in the vicinity of the droplet makes the hologram formation model different from the classically used diffraction pattern model of opaque spheres. These differences are mostly located in the low frequency part of the signal, in the central part of the diffraction pattern. We proposed to use an exclusion mask to reject the part of the signal that can not be correctly described with our diffraction pattern model. Using such an exclusion mask, we have shown that in-line digital holography employing the inverse problems approach has a precision that is sufficient to study the phenomenon of evaporation. The evaporation rate estimated from the digital holography measurements is consistent with the one found from a simple evaporation model. To reduce the processing time of the reconstruction method, we have taken advantage of the smooth movement of droplets in space and replaced the exhaustive search step of the algorithm by a parameter prediction step. This prediction step reduces the time costs of the iterative algorithm to 30 seconds per particle, improving the processing time by a factor of 10.

As mentioned above, a part of the droplet signature is caused by evaporation and is not explained by the diffraction model of non-evaporating droplets. A future work can be directed towards using a more rigorous light scattering model (e.g., generalized Lorenz-Mie theory for multi-layered sphere [31] ), to take the vapor film effect into account in the iterative algorithm. Such a model could provide additional information about the evaporation process.

Finally the prediction step could be further improved in the case of more complex trajectories using Kalman filters.

Appendix A

This appendix shows that, in the case of a spherical reference beam (produced by a point source located on the optical axis at distance z _s ), the mathematical model of hologram formation:

I_{holo} (x, y) = I_{ref}^{0} (x, y) - 2 ℜ {\underset{1}{\underset{︸}{A_{ref}^{0 *} (x, y)}} \cdot \underset{2}{\underset{︸}{A_{ref}^{z_{n}} (x_{n}, y_{n})}} \cdot \underset{3}{\underset{︸}{(ϑ_{n} * h_{z_{n}}) (x, y)}}}

can be approximated as

I_{holo} (x, y) \approx I_{ref}^{0} (x, y) + α_{n} \sqrt{I_{ref}^{0} (x, y)} . \frac{π r_{e}^{2}}{λ z_{e}} J_{1 c} (\frac{2 π r_{e} ρ_{e}}{λ z_{e}}) sin (\frac{π ρ_{e}^{2}}{λ z_{e}})

introducing a magnification factor

m = \frac{z_{s}}{z_{s} - z_{n}}

by means of a new parameter set θ _e = m θ _n with

ρ_{e} = \sqrt{{(x - x_{e})}^{2} + {(y - y_{e})}^{2}}

.

For simplicity, we consider the particle to be located on the optical axis i.e., x _e = y _e = 0 (see [36] for the general case). Let us consider, the terms 1, 2 and 3 of equation Eq. (20).

The first term corresponds to the reference wave amplitude $A_{ref}^{0 *} (x, y)$ on the hologram plane. The modulus of this term can be non-uniform. However it is equal to the square root of the background intensity: $I_{ref}^{0} (x, y)$ . The phase of this term is non-uniform too and depends on laser point source location. Under Fresnel approximation and considering only the phase modulation in x and y directions, the first term of Eq. (20) can be expressed as:

A_{ref}^{0 *} (x, y) = \sqrt{I_{ref}^{0} (x, y)} \cdot e^{- \frac{i π ρ_{e}^{2}}{λ z_{s}}} .

The second term is the amplitude of the reference wave across the particle’s aperture $A_{ref}^{z_{n}} (x_{n}, y_{n})$ . As the particle is located on the optical axis, the phase of this term is zero (we omit the phase modulation in z direction). We therefore obtain:

A_{ref}^{z_{n}} (x_{n}, y_{n}) = A_{ref}^{z_{n}} (0, 0) = A_{0}^{z_{n}} .

The third term is the object’s diffraction pattern. Assuming $r_{n} ≪ \sqrt{λ z_{n}}$ , this term can be written as (see [28 ] ):

[ϑ_{n} * h_{z_{n}}] (ρ) \approx \frac{π r_{n}^{2}}{i λ z_{n}} J_{1 c} (\frac{2 π r_{n} ρ_{e}}{λ z_{n}}) e^{\frac{i π ρ_{e}^{2}}{λ z_{n}}}

By using these simplifications in Eq. (20), we obtain:

I_{holo} (x, y) = I_{ref}^{0} (x, y) - 2 \sqrt{I_{ref}^{0} (x, y)} \cdot A_{0}^{z_{n}} \cdot ℜ {e^{- \frac{i π ρ_{e}^{2}}{λ z_{s}}} \cdot \frac{π r_{n}^{2}}{i λ z_{n}} J_{1 c} (\frac{2 π r_{n} ρ_{e}}{λ z_{n}}) \cdot e^{\frac{i π ρ_{e}^{2}}{λ z_{n}}}} .

The product of the two chirp functions can be simplified to $e^{- \frac{i π ρ_{e}^{2}}{λ z_{s}}} e^{\frac{i π ρ_{e}^{2}}{λ z_{n}}} = e^{\frac{i π ρ_{e}^{2}}{λ z_{e}}}$ , with z _e = mz _n. Defining $α_{e} = - \frac{2 A_{0}^{z_{n}}}{m}$ , we obtain

I_{holo} (x, y) = I_{ref}^{0} (x, y) + α_{e} \sqrt{I_{ref}^{0} (x, y)} \cdot ℜ {\frac{π r_{e}^{2}}{i λ z_{e}} J_{1 c} (\frac{2 π r_{e} ρ_{e}}{λ z_{e}}) e^{\frac{i π ρ_{e}^{2}}{λ z_{e}}}} .

Appendix B

In this appendix, the estimated range for the theoretical value of the evaporating rate K is detailed. When the evaporation law is formulated in terms of squared radius, the evaporation rate is given by

K = - (\frac{ρ_{g}}{ρ_{d}} D_{d} S_{h} ln (B_{M} + 1))

where ρ _g is the density of the gas around the droplet, ρ _d = 731.8 [ kg/m ³ ] is the density of of the diethyl ether, D _d = 0.918 × 10 ⁻⁵ [ m ² / s ] is the diffusion coefficient of the diethyl ether (vapor) in the air, S _h ≈ 2 is the Sherwood number and B _M is the Spalding mass number.

Since the real composition of the gas around the droplet is unknown, this gas is assimilated for simplicity to air, and its density ρ _g is assimilated to the density of dry air, that is not very different from that of humid air for our experimental conditions. A simple approximation of ρ _g is given by the law of ideal gas:

ρ_{g} = ρ_{air} = \frac{P_{atm}}{R_{d a} T_{r}}

where P _atm = 101325[ Pa ], R _da = 287[ J/Kg/K ] is the constant for dry air, and T _r is the temperature of reference which can be approximated as the ambient temperature ( T _r = 25° C ).

The Spalding mass number B _M can be approximated as

B_{M} = \frac{Y_{d} - Y_{\infty}}{1 - Y_{d}}

where Y _e represents the mass fraction of the vapor of diethyl ether at the droplet surface and Y _∞ ≈ 0 is the mass fraction of the vapor of diethyl ether at infinity (far from the droplet). Y _d is related to the saturation pressure of vapor at the interface which depends on the surface temperature:

Y_{d} = P_{{sat}_{d}} / (P_{{sat}_{d}} + (P_{atm} - P_{{sat}_{d}}) (M_{d a} + / M_{v_{d}}))

where P _{sat_d} depends on T, M _da = 28.965338 × 10 ⁻³ [ kg/mol ] is the molar mass of the dry air, and M _{v_d} = 74.12 × 10 ⁻³ [ kg/mol ] is the molar mass of the diethyl ether.

The correlation relating P _{sat_d} as a function of T has been taken from the Cheric Data Bank. It is written as:

P_{{sat}_{d}} = 1000 exp (- 12.43790 log (T) - \frac{6340.514}{T} + 95.14704 + 1.412198 \times 10^{- 5} T^{2})

where T ≈ 253° K to 258° K. The theoretical value of K is therefore ∈ [−7.5, −6.2] × 10 ⁻⁹ [ m ² /s ].

Acknowledgments

This work was funded by the MORIN project (3D Optical Measurements for Research and INdustry), which is supported by the french government through the “Agence Nationale de la Recherche” (ANR) and the “Programme Avenir Lyon-Saint-Etienne” (PAL-SE). The authors would like to thank Rolf Clackdoyle for his fruitful comments on the manuscript.

References and links

1. M. Sommerfeld and H.-H. Qiu, “Experimental studies of spray evaporation in turbulent flow,” International Journal of Heat and Fluid Flow 19, 10–22 (1998). [CrossRef]

2. S. M. Skippon and Y. Tagaki, “ILIDS measurements of the evaporation of fuel droplets during the intake and compression strokes in a firing lean burn engine,” Technical Report 960830, SAE International, Warrendale, PA (1996).

3. N. Damaschke, H. Nobach, and C. Tropea, “Optical limits of particle concentration for multi-dimensional particle sizing techniques in fluid mechanics,” Experiments in fluids 32, 143–152 (2002). [CrossRef]

4. D. Sugimoto, K. Zarogoulidis, T. Kawaguchi, K. Matsuura, Y. Hardalupas, A. Taylor, and K. Hishida, “Extension of the compressed interferometric particle sizing technique for three component velocity measurements,” in “13th international symposium on applications of laser techniques to fluid mechanicsLisbon, Portugal,”, 26–29 (2006).

5. F. Verpillat, F. Joud, P. Desbiolles, and M. Gross, “Dark-field digital holographic microscopy for 3D-tracking of gold nanoparticles,” Optics Express 19, 26044–26055 (2011). [CrossRef]

6. L. Dixon, F. C. Cheong, and D. G. Grier, “Holographic deconvolution microscopy for high-resolution particle tracking,” Optics Express 19, 16410–16417 (2011). [CrossRef] [PubMed]

7. X. Zhang, I. Khimji, U. A. Gurkan, H. Safaee, P. N. Catalano, H. O. Keles, E. Kayaalp, and U. Demirci, “Lensless imaging for simultaneous microfluidic sperm monitoring and sorting,” Lab on a Chip 11, 2535–2540 (2011). [CrossRef] [PubMed]

8. L. Repetto, E. Piano, and C. Pontiggia, “Lensless digital holographic microscope with light-emitting diode illumination,” Optics letters 29, 1132–1134 (2004). [CrossRef] [PubMed]

9. C. P. Allier, G. Hiernard, V. Poher, and J. M. Dinten, “Bacteria detection with thin wetting film lensless imaging,” Biomedical optics express 1, 762–770 (2010). [CrossRef]

10. J. R. Fienup, “Coherent lensless imaging,” in “Imaging Systems,” (2010). [CrossRef]

11. T. C. Poon, T. Yatagai, and W. Juptner, “Digital holography-coherent optics of the 21st century: introduction,” Applied Optics 45, 821 (2006). [CrossRef]

12. J. Coupland and J. Lobera, “Special issue : Optical tomography and digital holography,” Measurement Science and Technology 19, 070101 (2008). [CrossRef]

13. M. K. Kim, Y. Hayasaki, P. Picart, and J. Rosen, “Digital holography and 3D imaging: introduction to feature issue,” Appied Optics 52, DH1 (2013). [CrossRef]

14. J. Gire, L. Denis, C. Fournier, E. Thiébaut, F. Soulez, and C. Ducottet, “Digital holography of particles: benefits of the ‘inverse problem’ approach,” Measurement Science and Technology 19, 074005 (2008). [CrossRef]

15. C. Fournier, L. Denis, E. Thiebaut, T. Fournel, and M. Seifi, “Inverse problem approaches for digital hologram reconstruction,” in “Proceedings of SPIE ,”, 80430S (2011). [CrossRef]

16. J. Sheng, E. Malkiel, and J. Katz, “Digital holographic microscope for measuring three-dimensional particle distributions and motions,” Applied optics 45, 3893–3901 (2006). [CrossRef] [PubMed]

17. J. Lu, J. P. Fugal, H. Nordsiek, E. W. Saw, R. A Shaw, and W. Yang, “Lagrangian particle tracking in three dimensions via single-camera in-line digital holography,” New Journal of Physics 10, 125013 (2008). [CrossRef]

18. D. Nguyen, D. Honnery, and J. Soria, “Measuring evaporation of micro-fuel droplets using magnified DIH and DPIV,” Experiments in Fluids 1–11 (2010).

19. S. H. Lee, Y. Roichman, G. R. Yi, S. H. Kim, S. M. Yang, A. van Blaaderen, P. van Oostrum, and D. G. Grier, “Characterizing and tracking single colloidal particles with video holographic microscopy,” Optics Express 15, 18275–18282 (2007). [CrossRef] [PubMed]

20. F. C. Cheong, K. Xiao, D. J. Pine, and D. G. Grier, “Holographic characterization of individual colloidal spheres’ porosities,” Soft Matter 7, 6816–6819 (2011). [CrossRef]

21. D. Chareyron, J. L. Marie, C. Fournier, J. Gire, N. Grosjean, L. Denis, M. Lance, and L. Mees, “Testing an in-line digital holography “inverse method” for the lagrangian tracking of evaporating droplets in homogeneous nearly isotropic turbulence,” New Journal of Physics 14, 043039 (2012). [CrossRef]

22. J. Fung, K. E. Martin, R. W. Perry, D. M. Kaz, R. McGorty, and V. N. Manoharan, “Measuring translational, rotational, and vibrational dynamics in colloids with digital holographic microscopy,” Optics express 19, 8051–8065 (2011). [CrossRef] [PubMed]

23. F. Soulez, L. Denis, C. Fournier, É. Thiébaut, and C. Goepfert, “Inverse-problem approach for particle digital holography: accurate location based on local optimization,” Journal of the Optical Society of America . A 24, 1164–1171 (2007).

24. 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,” Journal of the Optical Society of America. A, Optics, Image Science, and Vision 24, 3708–3716 (2007). [CrossRef] [PubMed]

25. H. Royer, “An application of high-speed microholography: the metrology of fogs,” Nouvelle Revue d’ Optique 5, 87–93 (1974). [CrossRef]

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

27. G. B. Parrent and B. J. Thompson, “On the fraunhofer (far field) diffraction patterns of opaque and transparent objects with coherent background,” Journal of Modern Optics 11, 183–193 (1964).

28. M. Seifi, C. Fournier, L. Denis, D. Chareyron, and J.-L. Marié, “Three-dimensional reconstruction of particle holograms: a fast and accurate multiscale approach,” Journal of the Optical Society of America A 29, 1808–1817 (2012). [CrossRef]

29. M. Seifi, C. Fournier, and L. Denis, “HoloRec3D”. http://labh-curien.univ-st-etienne.fr/wiki-reconstruction/index.php/Main_Page, [Online; accessed 12-April-2013].

30. M. Seifi, C. Fournier, and L. Denis, “HoloRec3D : A free Matlab toolbox for digital holography,” Hal:ujm (2012).

31. L. Mees, N. Grosjean, D. Chareyron, J-L. Marie, M. Seifi, and C. Fournier, “Evaporating Droplet Hologram Simulation for Digital In-line Holography set-up with Divergent Beam,” Journal of Optical Society of America A 30, 2021–2028 (2013). [CrossRef]

32. M. Seifi, L. Denis, and C. Fournier, “Fast and accurate 3D object recognition directly from digital holograms,” Journal of the Optical Society of America A 30, 2216–2224 (2013). [CrossRef]

33. C. K. Law, T. Y. Xiong, and C. Wang, “Alcohol droplet vaporization in humid air,” International Journal of heat and mass transfer 30, 1435–1443 (1987). [CrossRef]

34. C. K. Law, “Recent advances in droplet vaporization and combustion,” Progress in energy and combustion science 8, 171–201 (1982). [CrossRef]

35. W. H. Press, T. A. Saul, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C, The Art of Scientific Computing (Cambridge University Press, 1992), cambridge university press ed.

36. C. S. Vikram and M. L. Billet, “Some salient features of in-line fraunhofer holography with divergent beams,” Optik 78, 80–83 (1988).

Accurate 3D tracking and size measurement of evaporating droplets using in-line digital holography and “inverse problems” reconstruction approach

Abstract

1. Introduction

2. Description of the experimental setup

3. Automated estimation of particles parameters

3.1. The parametric diffraction pattern model

3.2. Non-uniform background removal

3.3. Iterative parameter estimation algorithm

3.4. Pixel region exclusion mask

3.5. Rough parameter prediction

4. Experimental results

4.1. Calibration

4.2. Non-evaporating droplets

4.3. Evaporating droplets

5. Conclusion

Appendix A

Appendix B

Acknowledgments

References and links

Supplementary Material (2)

Cited By

Figures (8)

Equations (31)

Optics Express