Improvement of the size estimation of 3D tracked droplets using digital in-line holography with joint estimation reconstruction

Figure 1 shows the in-line holographic setup used for 3D tracking of free-falling droplets. This configuration is similar to that of [32]. The light emitted from a Nd-YVO₄ solid state laser ($\lambda =532$ nm Millennia IIs Spectra Physics) is low pass filtered and expanded using a f  =  100 mm focal-length lens associated with a 50 μm a diameter pinhole. The illumination beam is a spherical diverging wavefront impinging the water droplet jet under investigation. The jet consists of free-falling mono-dispersed water droplets generated by a piezo-electric injector (MJT-AT-01 MicroFab Technologies jetting device, 60 μm orifice diameter) located at 992 mm from the diaphragm and 472 mm from the imaging sensor. Interference between the field diffracted by the droplets and the reference beam (the part of the illumination beam that does not interact with the droplets) are recorded on a $1280\times 800$ pixel CMOS sensor (Phantom V611) with a pixel pitch of 20 μm and operating at 6200 frames per second. It should be noted that no collection optics is used for the holographic recording. The pixel fill factor, representing the ratio of the active surface of the pixel to the physical surface of the pixel, can be accounted for in our imaging model according to [33] and is $\kappa ={{\kappa}_{x}}\times {{\kappa}_{y}}=0.56$. The axial position of the water droplet is estimated roughly (based on the distance from the injector to the sensor) to be in the ${{z}_{\text{min}}}=0.46$ m to ${{z}_{\text{max}}}=0.48$ m range. To measure image magnification, the injection plane is considered as a reference plane. Magnification in the ith plane G_i was calibrated using a reticle and expressed as ${{G}_{i}}={{1.5.10}^{-3}}\Delta {{z}_{i}}+1.4735$ (see [35] for details about the calibration procedure). Here, $\Delta {{z}_{i}}$ is the algebraic distance between the injection and the ith planes (counted as positive in the light propagation direction). The z range $\left[{{z}_{\text{min}}};{{z}_{\text{max}}}\right]$ verifies $z\gg 4{{r}^{2}}/\lambda$, with r the droplet radius. Under this assumption, the droplets can be viewed as opaque spherical objects and the Huygens–Fresnel integral, which gives the intensity recorded on the imaging sensor, has an analytical solution [36]. The intensity ${{I}_{z}}(x,y)$ recorded for a spherical opaque particle at a distance of z from the imaging sensor is thus given by

$$\begin{eqnarray}&&{{I}_{z}}(x,y)\propto 1-\frac{1}{\lambda z}{{\mathcal{F}}_{\frac{x}{\lambda z},\frac{y}{\lambda z}}}\left\{\vartheta \left(\xi,\eta \right)\right\}\sin \left[\frac{\pi}{\lambda z}\left({{x}^{2}}+{{y}^{2}}\right)\right]\\ &&+{{\left[\frac{1}{\lambda z}{{\mathcal{F}}_{\frac{x}{\lambda z},\frac{y}{\lambda z}}}\left\{\vartheta \left(\xi,\eta \right)\right\}\right]}^{2}},\end{eqnarray} \tag{ 1 }$$

where ϑ is the aperture function of the spherical droplet defined as

$$\begin{eqnarray}&&\vartheta \left(\xi,\eta \right)=\left\{\begin{array}{*{35}{l}} 1 \\ 0 \end{array}\begin{array}{*{35}{l}} \text{if}\sqrt{{{\left(\xi -x\right)}^{2}}+{{\left(\eta -y\right)}^{2}}}\leqslant r, \\ \text{otherwise} \end{array},\right. \end{eqnarray} \tag{ 2 }$$

and $\left(\xi,\eta \right)$ are the coordinates in the aperture (droplet) plane. The Fourier transform of ϑ can be analytically derived and is given by

$$\begin{eqnarray}&&{{\mathcal{F}}_{\frac{x}{\lambda z},\frac{y}{\lambda z}}}\left\{\vartheta \left(\xi,\eta \right)\right\}=\pi {{r}^{2}}\left(\frac{\lambda z}{2\pi r\sqrt{{{x}^{2}}+{{y}^{2}}}}\right){{J}_{1}}\left(\frac{2\pi r\sqrt{{{x}^{2}}+{{y}^{2}}}}{\lambda z}\right),\end{eqnarray} \tag{ 3 }$$

where J₁ is the Bessel function of the first kind. It should be noted that the squared term in equation (1) is negligible in our experimental conditions as $\pi {{r}^{2}}/\left(\lambda z\right)\ll 1$. In the case of a diluted sample, optical interaction between particles can be disregarded. Therefore, multiple particles can be considered with an additive intensity model given by

$$\begin{eqnarray}&&{{I}_{z}}(x,y)\propto 1-\underset{i=1}{\overset{{{N}_{\text{p}}}}{\mathop \sum}}\,\frac{1}{\lambda {{z}_{i}}}{{\mathcal{F}}_{\frac{x}{\lambda z},\frac{y}{\lambda z}}}\left\{{{\vartheta}_{i}}\left(\xi,\eta \right)\right\}\sin \left[\frac{\pi}{\lambda {{z}_{i}}}\left({{x}^{2}}+{{y}^{2}}\right)\right],\end{eqnarray} \tag{ 4 }$$

where ${{N}_{\text{p}}}$ is the number of particles to be considered and i is the particle index. It should be noted that intrinsic in-line holographic experiment limitations apply to our reconstruction method. The shadow density has to be so that no more than 5–10% of the sensor area is occupied by the objects [34]. Moreover, the more the particle number, the less the hologram SNR and therefore the less the reconstruction accuracy.

Zoom In Zoom Out Reset image size

An example of a hologram obtained with the experimental setup shown in figure 1 is depicted in figure 2. As it can be seen, the contrast of the interference pattern is weak and the background is not spatially uniform. To avoid these problems, the acquired hologram sequences are pre-processed before being reconstructed.

Zoom In Zoom Out Reset image size

The aim of pre-processing the hologram is to enhance the signal to noise ratio (SNR) of the diffraction pattern. The principles of hologram pre-processing are illustrated in figure 3. Using the recorded hologram sequence, it is possible to build a background image by calculating the median of the sequence. Doing so makes it possible to obtain a background image with a noise structure much more representative of the hologram sequence than that of an 'empty' hologram recorded prior or after the image acquisition. For this purpose, the larger the number of acquired frames, the better the background. Here, all 1000 holograms were used to construct the background image. The constructed background image is then subtracted from the original hologram sequence, leading to an average removed sequence. It can be seen in figure 3 that the resulting sequence exhibits better fringe visibility and background uniformity. Moreover, the background noise (i.e. region of the hologram without interference signal) strongly resembles Gaussian white noise. This aspect will be very useful for the implementation of the proposed reconstruction scheme. This point is discussed in detail in section 4.

Zoom In Zoom Out Reset image size

Once the pre-processing step was complete, the hologram sequence is ready for hologram reconstruction. The proposed method is not based on backpropagation scheme as is often the case, but relies on inverse problem approaches. The next section is therefore devoted to the description of the reconstruction scheme.

Tracking is performed on a sequence of N  =  1000 holograms. In order to reduce the data processing time, values of z and r will be optimized in a range of  ±10% around roughly estimated values. More generally, the reconstruction approach remains valid as long as our imaging model remains valid (i.e. $z\gg 4\pi {{r}^{2}}/\lambda$). Seven to eight particles can be tracked on each hologram. For statistical reasons, only particle trajectories that span at least 200 holograms were considered. Moreover the holograms are cropped in order not to consider objects that exit from the injector (these are not spherical and do not corresponds to our image formation model). Therefore, according to the acquired holograms, ${{N}_{\text{traj}}}=15$ trajectories will be processed. For each trajectory, random stacks of ${{N}_{\text{h}}}=9$ holograms are built for the purpose of joint optimization.

An example of a hologram sequence is shown in figure 5(a) (see Media 1 for the dynamic sequence (stacks.iop.org/MST/27/045001/mmedia)). Here, white crosses corresponds to the position of the object detected using the first main steps of the reconstruction algorithm (i.e. classical IP reconstruction). A three-dimensional representation of the ${{N}_{\text{traj}}}=15$ is plotted in figure 5(b). The coordinates $(x,y)=(0,0)$ corresponds to the center of the hologram, while the injector is positioned at the top of the hologram (see figure 5(a)) and in the top left hand corner of figure 5(b). As it can be seen in figure 5(b), the extracted 3D trajectories are only weakly dispersed in both x and z directions. This can be explained by the fact that the imaging sensor is positioned to record the region of droplet injection, limiting the effect of the experimental environment on the droplet trajectories. To assess the accuracy of the 3D trajectory reconstruction, second order polynomial trajectories were fitted to the experimental data [43]. This choice was driven by the free-falling nature of the droplets. Results obtained for the ${{N}_{\text{traj}}}=15$ trajectories are given in table 1 in the 3rd to 5th columns for the classical IP, and in the 7th to 9th columns for joint reconstruction. From these results, it can be seen that the joint reconstruction process weakly influences the accuracy of the 3D tracking. This result was to be expected as neither x_i, nor y_i, nor z_i are constant throughout the image sequence. Moreover, it shows the low correlation between r and the $\left\{x,y,z\right\}$ parameters. It should be noted that the average processing time for each trajectory (classical IP approach and joint IP reconstruction for 8 particles tracked over 1000 holograms) is 8 h using Matlab 2015b on a high-end labtop computer (Windows 7-64 bit, Intel Core i7 3740QM @2.7 GHz with 24 Go RAM).

Zoom In Zoom Out Reset image size

Table 1. Data statistics for particle tracking and size estimation using both IP reconstruction and joint reconstruction.

Tracka	${{N}_{\text{hpt}}}~{{\left(\lfloor \frac{{{N}_{\text{hpt}}}}{{{N}_{h}}}\rfloor \right)}^{\text{b}}}$	$\sigma _{x}^{\text{IP}}~\left(\mu \text{m}\right)$	$\sigma _{y}^{\text{IP}}~\left(\mu \text{m}\right)$	$\sigma _{z}^{\text{IP}}~\left(\mu \text{m}\right)$	$\sigma _{r}^{\text{IP}}~\left(\mu \text{m}\right)$	$\sigma _{{\bar{r}}}^{\text{IP}}~{{\left(\text{nm}\right)}^{\text{c}}}$	$\sigma _{x}^{\text{Joint}}~\left(\mu \text{m}\right)$	$\sigma _{y}^{\text{Joint}}~\left(\mu \text{m}\right)$	$\sigma _{z}^{\text{Joint}}~\left(\mu \text{m}\right)$	$\sigma _{r}^{\text{Joint}}~\left(\mu \text{m}\right)$	$\sigma _{{\bar{r}}}^{\text{Joint}}~{{\left(\text{nm}\right)}^{\text{d}}}$
1	235 (26)	0.11	0.14	19.64	0.26	16.8	0.11	0.13	16.39	0.052	10.3
2	321 (35)	0.08	0.23	17.66	0.29	16.6	0.08	0.21	22.24	0.039	6.6
3	389 (43)	0.10	0.37	15.61	0.27	13.9	0.10	0.37	21.29	0.035	5.4
4	459 (51)	0.28	0.52	19.69	0.26	12.3	0.28	0.57	19.18	0.049	7.1
5	512 (56)	0.45	1.77	19.31	0.25	10.9	0.44	1.71	18.39	0.039	5.2
6	513 (57)	1.13	6.61	24.09	0.25	11.1	1.13	6.61	16.42	0.039	5.1
7	503 (55)	3.49	2.89	20.89	0.24	10.5	3.49	2.88	19.24	0.046	6.2
8	501 (55)	7.75	7.06	20.88	0.24	10.4	7.75	9.11	19.78	0.047	6.4
9	505 (56)	7.45	8.76	21.79	0.26	11.5	7.57	10.13	23.14	0.041	5.4
10	507 (56)	7.29	11.35	20.83	0.25	11.3	7.42	12.69	19.13	0.043	5.8
11	504 (56)	6.16	6.71	19.38	0.24	11.0	6.16	6.75	19.76	0.031	4.2
12	478 (53)	4.21	4.43	20.96	0.19	8.9	4.29	4.68	21.42	0.034	4.7
13	407 (45)	2.84	2.95	18.33	0.11	5.6	2.91	3.03	22.63	0.031	4.5
14	338 (37)	1.43	1.54	20.64	0.12	6.8	1.48	1.98	26.06	0.024	4.1
15	264 (29)	3.54	3.94	21.29	0.13	8.2	3.54	3.87	25.56	0.037	7.1

^aStandard deviations on the x, y, z particle positions ${{\sigma}_{x,y,z}}$ are estimated by fitting a theoretical free-falling droplet trajectory. ^b${{N}_{\text{hpt}}}$ corresponds to the number of hologram per trajectory. ^cEstimated over the number of holograms of the trajectory. ^dEstimated over the number of hologram stacks for each trajectory as given.

This processing time can be improved through parallel processing strategies. Nevertheless, these discussions fall out of the scope of this article. According to the values listed in table 1, the tracking precision (assumed to be the standard deviation between estimated positions and the fit of the trajectory, under the assumption of a second order polynomial trajectory [43]), is found to be ${{\sigma}_{x}}\times {{\sigma}_{y}}\times {{\sigma}_{z}}=3\times 3\times 20$ μm³, corresponding to $0.15\times 0.15\times 1$ pixels. It should be noted that we do not assess the absolute particle position (note that theoretical localization accuracy can be as low as one hundredth of a pixel [16]) but the precision of the trajectory fit.

As the radius is assumed to remain constant over the hologram sequence, its estimate should benefit from our joint reconstruction approach. The accuracy of our radius estimation will be assessed by two quantities: the standard deviation of the estimated radius distribution ${{\sigma}_{r}}$, and the standard error of the mean ${{\sigma}_{{\bar{r}}}}$. From the standard deviations of both IP and joint reconstruction estimation of the radius $\sigma _{r}^{\text{IP}}$, and $\sigma _{r}^{\text{Joint}}$, standard errors of the mean are defined as

$$\begin{eqnarray}&&\sigma _{{\bar{r}}}^{\text{IP}}=\frac{\sigma _{r}^{\text{IP}}}{\sqrt{{{N}_{\text{hpt}}}}},\end{eqnarray} \tag{ 5 }$$

for the standard error of the mean of r estimated by classical IP reconstruction. Here ${{N}_{\text{hpt}}}$ is the number of holograms per trajectory. The standard error of the mean (used to assess the accuracy of the estimation of the average) of the radius estimation using joint reconstruction is given by

$$\begin{eqnarray}&&\sigma _{{\bar{r}}}^{\text{Joint}}=\sigma _{r}^{\text{Joint}}\sqrt{\lfloor \frac{{{N}_{\text{h}}}}{{{N}_{\text{hpt}}}}\rfloor},\end{eqnarray} \tag{ 6 }$$

where $\lfloor \frac{{{N}_{\text{h}}}}{{{N}_{\text{hpt}}}}\rfloor$ is the integer part of the ratio of the number of holograms per stack to the number of holograms per trajectory.

The statistical results of the estimation of the radius using the two reconstruction schemes are listed in table 1. In addition, to better underline the interest of our joint reconstruction approach, changes in the radius estimation and statistical distribution of the droplet radius are illustrated in figure 6. These graphical representation were performed considering the 11th trajectory. In both figures 6(a) and (b), results obtained with classical IP reconstruction are plotted in blue, while results obtained with the joint reconstruction approach are plotted in orange. The interest of the random hologram stack construction is clear in figure 6(a). As can be seen, the estimated values of the radius are biased at the end of the hologram sequence and noise is clearly visible on the estimation. This behavior is no longer noticeable in the joint reconstruction results where the estimated values are barely scattered around the average estimate. This point is confirmed by the statistical distribution in figure 6(b). The results for all the trajectories are given in Media 2 (stacks.iop.org/MST/27/045001/mmedia). Note that estimates of the radius are less dispersed around the distribution average when a joint reconstruction scheme is used.

Zoom In Zoom Out Reset image size

Results obtained for the 15 tracks are summarized in table 1. From these 15 independent trajectories, and considering that the injector generate monodispersed droplet distributions, we can estimate the average standard deviations and standard errors for the two reconstruction methods. As the number of holograms differs from one trajectory to the next, weighted average values are calculated considering

$$\begin{equation} \langle\sigma_r^{\rm IP}\rangle=\sqrt{\frac{\sum_{i=1}^{15} N_{{\rm hpt}_i}{\sigma_r^{\rm IP}}^2_i}{\sum_{i=1}^{15} N_{{\rm hpt}_i}}} \end{equation} \tag{ 7 }$$

for the classical IP reconstruction, and

$$\begin{equation}\langle\sigma_r^{\rm Joint}\rangle=\sqrt{\frac{\sum_{i=1}^{15} \left\lfloor \frac{N_{{\rm hpt}_i}}{N_{\rm h}}\right\rfloor{\sigma_r^{\rm Joint}}^2_i}{\sum_{i=1}^{15} \left\lfloor \frac{N_{{\rm hpt}_i}}{N_{\rm h}}\right\rfloor}}\end{equation} \tag{ 8 }$$

for the joint reconstruction. Weighted averages on the standard errors were obtained the same way.

Therefore using classical IP reconstruction the radius can be estimated with a standard deviation $\langle \sigma _{r}^{\text{IP}}\rangle =0.235$ μm and a standard error $\langle \sigma _{{\bar{r}}}^{\text{IP}}\rangle =11.2$ nm, while using a joint reconstruction scheme these values are $\langle \sigma _{r}^{\text{Joint}}\rangle =0.0397$ μm and $\langle \sigma _{{\bar{r}}}^{\text{Joint}}\rangle =5.8$ nm resulting in a gain of 5.9 on the standard deviation and 1.9 on the standard error. It should be noted that in the case of white Gaussian noise, if the estimation of the radius computed using classical IP were averaged over ${{N}_{\text{h}}}=9$ holograms instead of performing the joint estimation reconstruction, the maximal gain in the radius estimation would be $\sqrt{{{N}_{\text{h}}}}=3$, which is half of the gain obtained using our method. Compared to a classical back-propagation reconstruction approach, the IP method has been shown to result in an improvement of the 3D localization of parametric objects [22]. As our approach is based on model-fitting, single point resolution has to be considered [16]. However, two-point resolution (e.g. Rayleigh criterion) is often considered.

For the purpose of comparison, and according to our experimental conditions, the Rayleigh criterion gives ${{\delta}_{x,y}}=44$ μm, and ${{\delta}_{z}}=3.6$ mm considering a numerical aperture of $\text{NA}=0.012$, which are 40 times higher for lateral resolution and 300 times higher in axial resolution.

Despite their consistency, the presented results can be biased by an inadequacy of the imaging model used for hologram reconstruction. To test the validity of our model, we performed particle hologram simulations with an electromagnetic model (Lorenz–Mie theory). Comparison with our simplified imaging model shows a difference of less than 2% in intensity. To estimate the bias thus introduced, simulations are realized under realistic conditions (parameters corresponding to our experiments, addition of a white Gaussian noise with the same variance as that of figure 5(a)). Reconstruction of the simulated holograms is realized using IP reconstruction with the Tyler/Thompson image formation model (see [36]). From the processing of these holograms, the bias introduced by the simplified model, is found to be about 1.5% on the z, r parameters estimation, while x, y estimation is not affected. This point shows the limits of our approach using a simplified model. Let us note that the whole process is 10 times longer using a Lorenz–Mie model due to its computation complexity. However this model can be used in a calibration step for the estimation of the reconstruction bias, making it possible to unbias measurements. These results were obtained in the particular case of digital in-line holography with parametric objects. It should however be noted that more complicated objects can be considered using a regularized MAP reconstruction approach. Moreover, the proposed method is not limited to digital in-line holography but can be generalized to any technique that is based on a reliable physical model.