# Quasi real-time analysis of mixed-phase clouds using interferometric out-of-focus imaging: development of an algorithm to assess liquid and ice water content

Pascal Lemaitre, M. Brunel, A. Rondeau, E. Porcheron, G. Gréhan

► To cite this version:

Pascal Lemaitre, M. Brunel, A. Rondeau, E. Porcheron, G. Gréhan. Quasi real-time analysis of mixed-phase clouds using interferometric out-of-focus imaging: development of an algorithm to assess liquid and ice water content. Measurement Science and Technology, 2015, 26 (12), pp.125403. 10.1088/0957-0233/26/12/125403. hal-01612365

HAL Id: hal-01612365

https://hal.science/hal-01612365

Submitted on 29 Sep 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

PAPER • OPEN ACCESS

Quasi real-time analysis of mixed-phase clouds
using interferometric out-of-focus imaging:
development of an algorithm to assess liquid and
ice water content

To cite this article: P Lemaitre *et al* 2015 *Meas. Sci. Technol.* **26** 125403

View the [article online](#) for updates and enhancements.

You may also like

- Defect data image enhancement method based on all-focus imaging algorithm
Yun Xie, Lujing Zhou, Xiaobin Zhang et al.

- Focus stacking in non-contact dermoscopy
Lennart Jütte, Zhiyao Yang, Gaurav
Sharma et al.

- Interferometric laser imaging for in-flight cloud droplet sizing
Christina Dunker, Christoph Roloff and
Arne Grassmann
# Quasi real-time analysis of mixed-phase clouds using interferometric out-of-focus imaging: development of an algorithm to assess liquid and ice water content

P Lemaitre¹, M Brunel², A Rondeau¹,², E Porcheron¹ and G Gréhan²

¹ Institut de Radioprotection et de Sûreté Nucléaire (IRSN), PSN-RES, SCA, LECEV, Gif-sur-Yvette 91192, France

² UMR 6614 CORIA, Université de Rouen, Avenue de l'université BP 12, 76801 Saint Etienne du Rouvray Cedex, France

E-mail: pascal.lemaitre@irsn.fr

Received 23 February 2015, revised 26 August 2015

Accepted for publication 28 August 2015

Published 26 October 2015

<figure><img src="image_3.png" /></figure>

CrossMark

## Abstract

According to changes in aircraft certifications rules, instrumentation has to be developed to alert the flight crews of potential icing conditions. The technique developed needs to measure in real time the amount of ice and liquid water encountered by the plane. Interferometric imaging offers an interesting solution: It is currently used to measure the size of regular droplets, and it can further measure the size of irregular particles from the analysis of their speckle-like out-of-focus images. However, conventional image processing needs to be speeded up to be compatible with the real-time detection of icing conditions. This article presents the development of an optimised algorithm to accelerate image processing. The algorithm proposed is based on the detection of each interferogram with the use of the gradient pair vector method. This method is shown to be 13 times faster than the conventional Hough transform. The algorithm is validated on synthetic images of mixed phase clouds, and finally tested and validated in laboratory conditions. This algorithm should have important applications in the size measurement of droplets and ice particles for aircraft safety, cloud microphysics investigation, and more generally in the real-time analysis of triphasic flows using interferometric particle imaging.

Keywords: interferometric out-of-focus imaging, drop characterisation, ice crystal characterisation, freezing conditions

(Some figures may appear in colour only in the online journal)

## Introduction

Since the 1990s, aircraft engine incidents and abnormalities in aerodynamic speed and temperature measurements have occurred at high altitude and low temperature. They were attributed to the icing of both the engines and the speed and

<figure><img src="image_4.png" /></figure>

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any

further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

temperature sensors. The detailed analysis of the parameters recorded during these events showed that the temperatures recorded were mostly lower than -40 °C and no ice formation was detected on the aircraft structure.

These findings led to the idea that a so far unidentified phenomenon occurred, different from conventional icing caused by ice accretion resulting from the impact of supercooled droplets. This kind of conventional icing is already taken into account in the aircraft certification rules. However, it appears that the ingestion of large amounts of ice micro-crystals at
high altitude is the cause of both engine malfunction and speed indication abnormalities. Changes in the aircraft certification rules are thus aimed at considering these risks related to micro-crystals. It is necessary to develop on-board instruments for the detection of air masses that may be hazardous due to the simultaneous presence of micro-crystals and supercooled droplets. Different airborne instruments such as the Forward Scattering Spectrometer Probe (FSSP, De Araujo Coelho et al 2005) and Phase Doppler Interferometry (PDI, Chuang et al 2008) can provide the measurement of the droplet size distribution in clouds. However, they do not give any information on the ice crystal sizes. Only the Cloud Particle Imager (CPI, Lawson et al 2006) and Small Ice Detector 3 (SID-3, Ulanowski et al 2014) can provide the size distributions of both droplets and ice crystals. However, for the CPI probe the measurements are not continuous and there is no automatic processing of the images. The SID-3 can only measure 30 particles per second. That seems an insufficient sample to alert in case of icing conditions.

The technique selected to achieve this goal is interferometric out-of-focus imaging. This technique, first introduced by König et al (1986) and then improved by Glover et al (1995), is conventionally used in laboratories for droplet-size measurements. An optical sensor using this technique was recently developed by Porcheron et al (2015). It provides information on the droplets probability density function in warm clouds for a better understanding of their formation mechanisms. However, the image processing of this probe (Quérel et al 2010) fails in the presence of ice crystals and does not give any information on the droplet concentration. Recent works, carried out by Brunel et al (2014a, 2014b), showed that this technique could be extended to characterise irregular particles, such as ice crystals.

Based on these studies, we seek to develop a new airborne instrument to measure in real time the concentrations per size class of both droplets and ice crystals encountered by the aircraft. The final objective is to determine, in real time, the liquid and ice water contents (LWC and IWC). These two parameters respectively characterise the amounts of liquid water and ice in suspension in the cloud (generally expressed in g m⁻³). Planes' certification rules are based on these two parameters. This article focuses on the image analysis method selected to determine the IWC and LWC. The processing of the images should be at a rate compatible with the in-flight detection of icing conditions and thus exceed 10 Hz. The processing of each image should consequently be shorter than 0.1 s.

# 1. General background

Interferometric out-of-focus imaging is a technique commonly used in laboratories to measure the size of spherical particles, droplets (Maeda et al 2000; Lemaitre et al 2007) or bubbles (Dehaeck and van Beeck 2007, Shen 2014). For a droplet, it is based on the analysis of the interference between the light rays reflected and those refracted by a transparent spherical particle. Indeed, when a droplet is illuminated by a coherent wave, a fraction of it will be reflected by the droplet

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Principle of the interferometric out-of-focus imaging technique.</figcaption></figure>

and another will be refracted. If this droplet is observed in the focal plane of an optical system, two ‘glare points’ (Van de Hulst and Wang 1991) associated with the reflected and refracted rays are respectively observed (figure 1). If a sensor is now placed outside the focal plane, interference fringes appear. The Lorenz–Mie theory and geometrical optics both enable the theoretical relationship between the frequency of the fringes and the size of the droplets to be determined.

Mounaïm-Rousselle and Pajot (1999) established a simple relationship that enables the frequency of the interference fringes ($F_{max}$) to be linked to the diameter of the droplet ($d$) at the origin of the interferogram (equation (1)).

$$d = 2\lambda F_{max} \left( \cos\left(\frac{\theta}{2}\right) + \frac{m \sin\left(\frac{\theta}{2}\right)}{\sqrt{m^2 - 2m \cos\left(\frac{\theta}{2}\right) + 1}} \right)^{-1} \quad (1)$$

In this equation, $\theta$ is the scattering angle (figure 1), m is the refraction index, and $\lambda$ is the wavelength of the incident wave. This relationship was validated experimentally by Lemaitre et al (2007) on trains of monodisperse water droplets.

More recently, Brunel and Shen (2013) developed a simplified model that enables interferometric out-of focus image formation to be simulated through an optical system. This approach is based on solving the generalised Huygens–Fresnel integral by considering the two ‘glare points’ in the source plane as two point emitters. Their positions, phase shifts, and relative intensities are calculated using geometric considerations and a Debye series expansion of the Lorenz–Mie theory (Debye 1908). The two interferograms presented in figure 1 outside the focal plane are calculated using this approach, by assuming that the optical system of collection consists of a thin lens with a focal distance ($f$) of 100 mm and a diameter ($\phi$) of 40 mm. These two interferograms have out-of-focus distances of, respectively, 250 and 600 µm. The image of the ‘glare points’ in the focal plane was obtained from Lemaitre et al (2007). This out-of-focus image simulation model was validated for both bubbles and droplets (Shen 2014).
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Adaptation of the interferometric out-of-focus imaging technique for ice crystals.</figcaption></figure>

The detection and characterisation of irregular rough particles, such as ice crystals, is much more complex. First of all because there is currently no rigorous theory that enables the scattering of light by such particles to be physically described, but also because there are many varieties of ice crystal sizes and shapes (Pruppacher and Klett 1998).

Ulanowski *et al* (2012) were the first to investigate the bi-dimensional structure of interferograms scattered by ice crystals (and other irregular particles). They showed the speckle structure of the forward scattering pattern and developed an instrument (Small Ice Detector-3, Ulanowski *et al* 2014), based on the bi-dimensional analysis of these interferograms. The approach used to understand the signals collected is based on the Fraunhofer approximation, assuming an arbitrary number of point emitters within the particle.

Brunel *et al* (2014a, 2014b) adapted their out-of-focus image simulation model based on the same hypotheses as those used by Ulanowski *et al* (2012). Thus, they modelled the ice crystal in the object plane by a set of N point emitters (Dirac emitters) with random phases, distributed arbitrarily within the geometrical boundaries of the ice crystal (figure 2). Then, solving the generalised Huygens–Fresnel integral for these N point emitters, it is possible to simulate the interferogram collected by an out-of-focus sensor, through modelling of the optical collection system using the transfer-matrix formalism.

Figure 2 shows the general outline of the method. In this figure the interferogram presented is directly simulated using the model by Brunel *et al* (2014a, 2014b); a speckle-like pattern is observed. Brunel *et al* (2014b) proposed measuring the speckle grain size by calculating the 2D autocorrelation function (Goodman 2009) of the interferogram. They thus showed that the width of the central peak of the 2D autocorrelation function is inversely proportional to the dimensions of the crystal. This approach was very recently validated experimentally by Brunel *et al* (2015). They further demonstrated that the 2D-autocorrelation of the particle shape is quantitatively

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Optical device for interferometric out-of-focus imaging.</figcaption></figure>

given by the 2D Fourier transform of the speckle pattern it scatters. This is performed by analyzing the speckle patterns scattered by irregular particles and measuring at the same time their shape by optical microscopy.

Thus, interferometric out-of-focus imaging allows the measurement of the size of ice crystals and water droplets over large fields. Moreover, an interferogram simulation formalism validated experimentally for these two types of particles is presently available.

## 2. Interferometric out-of-focus imaging setup for mixed phase characterisation

The optical setup for interferometric out-of-focus imaging is similar to a PIV device (Adrian and Yao 1985), apart from the sensor which is placed outside the focal plane. An optical assembly was developed in our laboratory and is presented in figure 3.

This optical assembly is very straightforward; a laser sheet is formed using a frequency-doubled Nd:YAG laser and a cylindrical lens. A camera is set to view this laser sheet at an angle ($\theta$) of 90°. The sensor used for this assembly is an sCMOS sensor (2560 × 2160 pixels², 16 bits per pixel). It is equipped with a Zeiss objective lens with a focal distance of 25 mm and a f/2.8 aperture. This lens is placed at a distance ($z1$) of 100 mm from the volume of measurement and the diaphragm is opened to the maximum, that is, an aperture angle ($\alpha$) of 5.1°. The field of view of this optical collection system is 10 cm × 10 cm. The laser sheet thickness is 1 mm. As a consequence, the volume of measurement is 10 cm × 10 cm × 1 mm.

This optical system can be modelled with the use of transfer matrices formalism. $M_1$ and $M_2$ are, respectively, the transfer matrices modelling the light path propagation in free space over the distances $z1$ and $z2 + \Delta p$ (figure 3). $M_L$ is the transfer matrix modelling the refraction by the lens. The total transfer matrix ($M_{\text{tot}}$) describing the propagation of the light path between the volume of measurement and the image sensor is the product of these three matrices (equation (2)).
<figure><img src="image_8.png" /><figcaption>Figure 4. Interferometric out-of-focus images of a binary mixture of droplets and ice crystals.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 5. Left image: Speckle interferogram. Middle image: Fringe interferogram. Right image: Preconditioned image for both speckle and fringe interferograms; $\vec{V}$ are gradient vectors.</figcaption></figure>

$$ M_{tot} = \begin{pmatrix} A_{tot} & B_{tot} \\ C_{tot} & D_{tot} \end{pmatrix} = M_2 \times M_L \times M_1 \\ = \begin{pmatrix} 1 & z1 \\ 0 & 1 \end{pmatrix} \times \begin{pmatrix} 1 & 0 \\ \frac{-1}{F} & 1 \end{pmatrix} \times \begin{pmatrix} 1 & z2 + \Delta p \\ 0 & 1 \end{pmatrix} \quad (2) $$

Figure 4 is an example of a set of images obtained with this assembly for a mixed phase. The droplets are produced using a nozzle and the crystals are produced by re-suspending ice initially accreted around a pipe cooled to -40 °C. In these images, each interferogram consisting of vertical fringes is a signal scattered by a droplet and the speckle signal is scattered by a cluster of ice crystals.

In order to deduce the LWC and IWC for a known volume of measurement, it is necessary to determine the respective sizes of the ice crystals and the droplets within this volume, at a rate compatible with the in-flight detection of icing conditions.

### 3. Real-time analysis of mixed clouds (droplets and ice crystals)

With the image global processing strategy developed by Quérel et al (2010), the LWC and IWC are not readily obtained. As far as it makes no distinction between the interferograms scattered by droplets and those scattered by ice crystals, it totally

<figure><img src="image_10.png" /><figcaption>Figure 6. Determination of the potential centres (+) and the cluster centre (×).</figcaption></figure>

fails in presence of mixed phase. It is necessary to develop a new algorithm to detect each interferogram in the image and to distinguish whether it was scattered by a droplet (vertical fringes, figure 1) or by an ice crystal (speckle, figure 2). Eventually, these signals are analysed using the appropriate method.

#### 3.1. Detection of each interferogram in the image

Detecting each interferogram (interference disk) in an out-of-focus image consists in identifying each disk in the image
<figure><img src="image_11.png" /><figcaption>Figure 7. Example of interferogram detection in the set of experimental images presented in figure 4.</figcaption></figure>

(centre and diameter). Two approaches are envisaged for this purpose: The *gradient pair vector* (GPV) method and the *circular hough transform* (CHT). The latter is much more conventional but very time-consuming in terms of calculation time, because it requires the construction of an accumulator matrix (Kimme et al 1975).

The vector pair method, introduced by Rad et al (2003) is, as its name indicates, based on the properties of the vectors normal to a disk. The first step consists in pre-conditioning the image in such a way that the inside of each disk is darker than the outside. Then, the gradient of the grey levels of the image is computed for each pixel. Then, the magnitude of this gradient is thresholded and finally the grey levels are inverted. Figure 5(γ) is the result of this pre-conditioning applied to, respectively, the speckle (figure 5(α) and fringe (figure 5(β)) interferograms.

Then, the gradient of this pre-conditioned image is calculated for each pixel (vectors in figure 5(γ)). For symmetry reasons, for each vector normal to a disk, there is another vector of the same magnitude and opposite direction ($\vec{V}_j$ and $\vec{V}_k$ in figure 5(γ)). These vectors are called vector pairs.

The second step of the algorithm by Rad et al (2003) consists in seeking all of the vector pairs of the image; these vectors must meet the following criteria:

$$ \left\{ \begin{array}{l} \frac{\overrightarrow{A_k A_j} \cdot \overrightarrow{V_j}}{\| \overrightarrow{A_k A_j} \| \cdot \| \overrightarrow{V_j} \|} \approx 1 \\ \\ \frac{\overrightarrow{V_k} \cdot \overrightarrow{V_j}}{\| \overrightarrow{V_k} \| \cdot \| \overrightarrow{V_j} \|} \approx -1 \end{array} \right. \qquad (3) $$

The different vectors of this system of the equation are defined in figure 5(γ). In our particular case, the radius of the interferograms is known (R), because it is related to the size of the lens and to the out-of-focus distance. Therefore, a third criterion based on the distance between the pairs of vectors is added:

$$ \overrightarrow{A_k A_j} \approx 2R \qquad (4) $$

Therefore, for each gradient vector $\vec{V}_j$ of the image, a vector $\vec{V}_k$ located at $A_k$ that meets the criteria of equations (3) and (4) is sought. If there is one such vector, a potential interferogram centre identified at the coordinate (Ox, Oy) corresponds to the middle of the segment $[A_k A_j]$ (+in figure 6).

Finally, a matrix of the triplets (Ox, Oy, R) is constructed and clustered using a density-based algorithm (Ester et al 1996, Daszykowski et al 2001) for identified potential centres. This finally provides the centre of each cluster (×in figure 6).

Applied to the set of images shown in figure 4, the results presented hereinafter are obtained (figure 7). On these images each circle corresponds to an interferogram detected with the GPV method.

In these few tests, the quality of the interferograms detection can be observed; even those with high overlap rates and very low contrast are perfectly identified (i.e. disk 2 in figure 7(β)).

The second step of the image processing consists in differentiating the speckle signals scattered by the ice crystals from the vertical fringes scattered by the droplets.

## 3.2. Differentiation of the speckle signals (scattered by ice crystals) from the vertical fringes (scattered by droplets)

This step would be very simple and very fast if all fringes of the interferograms scattered by the droplets were perfectly vertical. In that case the standard deviation of the grey levels in the vertical direction of the interferograms is at least 100 times smaller in the case of fringes compared to speckle. Unfortunately, as shown in figure 7 (interferograms β-4), the fringes are sometimes slightly slanted.

Brunel et al (2015) showed that the structure of the 2D Fourier transform of speckle patterns is correlated to the shape of the particle. Thus, a criterion is sought in the k-space that will enable the differentiation of speckle and fringe patterns. To do so, each overlapping zone of the interferograms is set to zero, and the 2D Fourier transforms of the resulting pattern is calculated. Figure 8 shows this result applied to Image α (presented in figure 7). The top line corresponds to the non-overlapping areas of each interferogram, and the bottom line presents the 2D Fourier transforms of those interferograms.
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Top line: non-overlapping areas of each interferogram of Image α; bottom line: 2D Fourier transform of each interferogram.</figcaption></figure>

In this figure, it appears that the number of spots in the k-space is directly correlated with the nature of the interferogram. The 2D Fourier transform of fringes has 3 spots close to the horizontal axis, while the speckle has a cloud of spots. The number of spots in the Fourier space thus tells the nature of the particle. This step is not time consuming because the 2D Fourier transform is needed in the next step of the algorithm (section 3.3) dedicated to determining the size of the droplets and ice crystals.

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Determination of the interference fringe frequency. In the image on the left, the parts overlapped by several interferograms are cropped. The figure on the right shows the associated 2D FFT.</figcaption></figure>

### 3.3. Determination of the size of the droplets and ice crystals

Once the speckle signals (scattered by ice crystals) are differentiated from the interference fringes scattered by droplets, each interferogram is processed appropriately.

#### 3.3.1. Analysis of the interferograms scattered by droplets.
For interferograms scattered by droplets, their diameter is calculated from equation (1). In that equation, the frequency of the fringes ($F_{max}$) is calculated using a 2D Fast Fourier Transform (calculated using the algorithm of Cooley and Tukey (1965)) on the areas of the interferogram with no overlapping. Figure 9 illustrates an example of the calculation of the fringe frequency. The example illustrated (interferogram No.2 in figure 7(β) is one of the most delicate configurations, because the contrast of the fringes is low and the interferogram is overlapped by two other figures. Despite this, all of the spots in the 2D FFT are easily identified.

#### 3.3.2. Analysis of the interferogram scattered by ice crystals.
For interferograms scattered by irregular particles like ice crystals, Brunel et al (2014b) suggested a simple relationship linking the speckle grain size ($\delta$) to the dimension of a particle ($\Delta$).

$$ \Delta = \frac{\lambda B_{tot}}{\delta} \qquad (5) $$

In this equation, the coefficient $B_{tot}$ is the coefficient B of the total transfer-matrix between the plane of the scattering particle and the sensor plane (equation (2)). The average speckle grain size ($\delta$) can be calculated using the width of the

central peak of the 2D autocorrelation of the speckle pattern (Goodman 2009).

In order to optimise the computation time, this 2D autocorrelation function ($R(x, y)$) is calculated using the Wiener-Khinchin theorem (equation (6); Wiener 1930).

$$ S(\nu) = \int_{-\infty}^{+\infty} R(x)e^{-i2\pi\nu x}dx = \hat{f}\{R(x)\} \qquad (6) \\ \text{With } R(X) = \hat{f}^{-1}\{S(f)\} $$

$S(f)$ is the power spectral density and $\hat{f}$ is the Fourier transform. This theorem is then applied in two dimensions. For illustration, this algorithm is used to calculate the autocorrelation of two speckles detected in images (α) and (γ). Finally, the speckle grain size $\delta$ is calculated by measuring the width of the autocorrelation peak at 70% of the maximum correlation (Brunel et al 2014b). The speckle grain is outlined in red in figure 10.

It is observed that the speckle grain is not circular. It is not surprising as the particles that scatter these signals are not circular either. Besides, the major and minor axis of the speckle grain can be determined. They are inversely proportional to the major axis and minor axis of the ice crystal (equation (7)).

$$ \Delta_{\min} = \frac{\lambda B_{tot}}{\delta_{\max}} \quad \text{and} \quad \Delta_{\max} = \frac{\lambda B_{tot}}{\delta_{\min}} \qquad (7) $$

Likewise, it is simple to measure the orientation of the ice crystal by measuring the angle between its major axis and
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Calculation of the autocorrelation functions of speckle interferograms: On the left autocorrelation of the speckle in disk 1 from figure 7(α) and on the right autocorrelation of the speckle in disk 1 from figure 7(γ).</figcaption></figure>

the horizontal axis. To assess the IWC, the average speckle grain size ($\bar{\delta}$) representative of the mass of the ice crystal is needed. This is performed by calculating the diameter of a disk of same surface as the speckle grain.

$$ \bar{\delta} = \sqrt{\frac{4 S_{\text{grain}}}{\pi}} \qquad (8) $$

In equation (8) $S_{\text{grain}}$ is the surface of the speckle grain. Then, the average diameter of the crystal ($\bar{\Delta}$) is calculated using the equation below.

$$ \bar{\Delta} = \frac{\lambda B_{\text{tot}}}{\bar{\delta}} \qquad (9) $$

Thus, knowing the dimensions of the volume of measurement of each image ($V = 10 \text{ cm} \times 10 \text{ cm} \times 1 \text{ mm}$) and the size of the droplets ($d$) and ice crystals ($\Delta$), the IWC and the LWC are determined (equation (10)).

$$ \begin{aligned} \text{LWC} &= \sum_{i=1}^{N_{\text{fringe}}} \frac{\pi \cdot d_i^3 \rho_{\text{drop}}}{6 \cdot V \cdot N_{\text{image}}} \\ \text{IWC} &= \sum_{i=1}^{N_{\text{speckle}}} \frac{\pi \cdot \bar{\Delta}_i^3 \rho_{\text{ice}}}{6 \cdot V \cdot N_{\text{image}}} \end{aligned} \qquad (10) $$

In equation (10) $N_{\text{image}}$ is the number of images on which the measurement is performed. $N_{\text{fringe}}$ and $N_{\text{speckle}}$ are, respectively, the numbers of interferograms detected and attributed by the algorithm to be scattered by droplets or ice crystals.

# 4. Validation and performances

## 4.1. Validation of the algorithm

The algorithm validation is carried out with synthetic images. Each synthetic image is composed on average of 10 interferograms. For each interferogram, the input data such as particle phase (droplet or ice crystal) and size are randomly selected (the Gaussian size distributions presented in table 1). These input data are then recorded, for further comparison with the algorithm results (interferogram by interferogram).

<table><caption>**Table 1.** Comparisons between the image data and the results of the algorithm.</caption><thead><tr><th></th><th>In the images generated</th><th>Algorithm result</th></tr></thead><tbody><tr><td>Total number of droplets</td><td>4984</td><td>4944</td></tr><tr><td>Total number of crystals</td><td>5016</td><td>5056</td></tr><tr><td>Average diameter of the droplet size distribution</td><td>100 µm</td><td>102 µm</td></tr><tr><td>Standard deviation of the droplet size distribution</td><td>50 µm</td><td>55 µm</td></tr><tr><td>Average diameter of the crystal size distribution</td><td>60 µm</td><td>64 µm</td></tr><tr><td>Standard deviation of the crystal size distribution</td><td>25 µm</td><td>30 µm</td></tr></tbody></table>

For interferograms scattered by ice crystals (detailed in section 1), a final random selection is performed in order to determine the number (between 14 and 30 per crystal) and location of the point emitters at the origin of the speckle. Finally, the interferogram scattered by each particle is calculated according to the generalised Huygens–Fresnel integrals associated with transfer matrices formalism to describe the imaging system (Brunel *et al* 2014a, 2014b).

The imaging system considered is simplified (relative to the one used in section 2). It consists of a thin lens with a focal distance ($F$) of 24 mm and a diameter of 8.5 mm (in order to have an aperture similar to the one in our experiments, $F/2.8 = 35.7$ mm). It is located at 90 mm from the volume of measurement ($z1 = 90$ mm). The incident wavelength is 532 nm. The virtual sensor is located 10 mm beyond the focal plane ($\Delta p = 10$ mm).

Once the interferograms have been computed, their locations are randomly selected and then they are projected onto a virtual sensor with a resolution of 2000 pixels × 2000 pixels and an 8-bit dynamic range.

Figure 11 shows an example of a synthetic image. Each circle in figure 11 corresponds to an interferogram detected with the GPV method. Thousands of images are thus synthesised numerically and processed with both GPV and CHT methods.
<figure><img src="image_15.png" /><figcaption>Figure 11. Example of a synthetic image processed with the GPV method.</figcaption></figure>

The results of both methods are very similar (the same detection ratio and interferogram location). However, the CHT is 13 times slower than the GPV.

Then, the nature of each interferogram is established (speckle or fringe signal, section 3.2). Finally the diameter of the particles is calculated (equation (1) for droplets and 9 for ice crystals). To compute the ice crystals' average diameter (Δ) the Btot coefficient of the total transfer matrix is needed. Btot is deduced from equation (2):

$$B_{\text{tot}} = z1 \left( 1 - \frac{(z2 + \Delta p)}{F} \right) + z2 + \Delta p = -0.0282 \quad (11)$$

And z2 is computed with the conjugation relation:

$$\frac{1}{z1} + \frac{1}{z2} = \frac{1}{F} \quad (12)$$

Table 1 and figure 12 compare the results of the random selections from thousands of synthetic images and the final result of the algorithm.

These results are globally very satisfactory; the interferogram detection rate is very high (over 99%). Only the interferograms in the corners of the images are undetectable. Indeed, there is no symmetrical vector in this configuration. The size distributions of the ice crystals and droplets are close to the ones imposed on the synthetic images.

A slight increase in the standard deviations of the particle size distributions is noted, compared to the input data of the synthetic images (around 10% for both ice crystals and droplets). This is also related to the fact that for some interferograms with high overlapping, such as interferograms 2 and 3 in figure 11, the area of each pattern that is not overlapped is smaller. However, the frequency resolution of 2D FFT is proportional to the number of pixels on which it is performed. As a consequence, for strongly overlapped interferograms, a loss of accuracy regarding the frequency is induced. Therefore, the precision on the droplet size measurement is also reduced.

It is also noticed that for 40 of the 4984 interferogram produced (0.8%), the algorithm analyses the fringes patterns (scattered by droplets) as speckle interferograms (scattered by ice crystals). This is observed on images when an

interferogram scattered by very small droplets (with small fringe frequency) is overlapped by another interferogram. For this particular case, the non-overlapping area is so small that it does not contain two fringes any more. Thus, when the differentiation between the ice crystals and the droplets is performed, the three spots do not appear clearly in the k-space and thus the particle is interpreted as a small ice crystal. As it is rarely observed, and only for the finest droplets, it has no measurable effect neither on the LWC and IWC deduced.

Finally, a slight increase in the average diameter of the droplets (10%) and crystals (6%) is also measured; this is observed in extreme configurations with major interferogram overlapping, as well as the non-overlapped areas of the interferograms being smaller than the fringe spacing (of the little droplets) and the speckle grain size (of the little crystals).

From the diameter, nature, and number of particles detected on each image (table 1 and figure 12) the LWCalgo and IWCalgo are calculated (equation (13)) and then compared to LWCsynim and IWCsynim computed from the input data of same synthetic images.

$$\begin{align}
\text{LWC}_{\text{algo}} &= \sum_{i=1}^{N_{\text{fringe}}} \frac{\pi \cdot d_i^3 \rho_{\text{drop}}}{6 \cdot V \cdot N_{\text{image}}} = 1.09 \text{ LWC}_{\text{synim}} \tag{13} \\
\text{IWC}_{\text{algo}} &= \sum_{i=1}^{N_{\text{speckle}}} \frac{\pi \cdot \bar{\Delta}_i^3 \rho_{\text{ice}}}{6 \cdot V \cdot N_{\text{image}}} = 1.4 \text{ LWC}_{\text{synim}} \nonumber
\end{align}$$

In these equations $d_i$ and $\bar{\Delta}_i$ are, respectively, the diameters of the droplets and ice crystals determined from the interferogram analysis, $N_{image}$ is the number of images on which the comparison is performed ($N_{image} = 1000$), $N_{fringe}$ and $N_{speckle}$ are, respectively, the numbers of interferograms detected and attributed by the algorithm to be scattered by droplets or ice crystals ($N_{fringe} = 4944$, $N_{speckle} = 5056$; table 1), and V is the volume of measurement of each image (10 cm × 10 cm × 1 mm).

The comparisons are in satisfactory agreement with a 10% overestimation of the LWC and 40% overestimation of the IWC. Note that for aircraft security reasons it is crucial to not underestimate the LWC and IWC.

As the main source of error identified is overlapping between the interferograms, the same comparisons are performed while increasing the number of interferograms per image from 5 to 25, with the same optical configuration (i.e. the same diameter of the interferograms) and the same initial size distributions for both ice crystals and droplets. The result of this comparison is presented in figure 13.

An increase in both the IWC and LWC overestimation is progressively observed up to a concentration of 20 interferograms per image. This is due to the growth of the overlapping ratio of the interferograms leading to a loss of precision on the determination of particle diameter. Over 20 interferograms per images, the GPV method does not detect all the interferograms any more due to increasing overlapping. This induces a decrease in both the IWC and LWC.

Such computations will help to optimise the final instrument design; especially the degree of defocusing that induces
<figure><img src="image_16.png" /><figcaption>**Figure 12.** Comparison of the size distributions on the synthetic images with the results of the algorithm for both droplets and ice crystals.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 13.** Influence of the number of interferograms in the image on the IWC and LWC precision.</figcaption></figure>

the size of the interferograms and thus drives their degree of overlapping. Indeed, regarding icing certification rules, it will be important to minimise the interferogram overlapping for such particle concentrations. Moreover, the pixel size of the camera is another important parameter. Indeed, for the same interferogram dimensions, decreasing the pixel size induces an increase in the number of pixels on which the FFT is performed. As a consequence, the final precision of the particle diameter measurement (for both droplets and ice crystals) is improved.

Finally, the algorithm is tested on a sample of 100 experimental images (of which three examples are presented in figure 4), with the experimental device described in section 2. The overall result of the interferogram detection is similar to the one presented in figure 7, and is thus very satisfactory. Contrary to the synthetic images analysis, no inversion in the nature of the particle is identified.

Unfortunately, it was not possible to perform further experimental validations. Indeed, no other instrument is available in our laboratory with the capacity (like interferometric out-of-focus imaging) to measure the size of ice crystals.

## 4.2. Performances of the algorithm

The complete algorithm was coded in Matlab on an office laptop (Intel Pentium, i5, 1.17 GHz). The average calculation time is around 0.15 s/image for interferogram detection using the GPV method (described in section 3.1), while the CHT detection method, which is much more conventional, takes around 2 s/image. Then, each interferogram analysed is around 0.05 s. Thus, the analysis time for an image depends on the number of interferograms per image.

## Conclusion

A method for the real-time analysis of images produced using the interferometric out-of-focus imaging technique has been developed. Thanks to the recent studies by Brunel et al (2014a, 2014b, 2014c), this technique was extended to the analysis of interferograms obtained from rough particles, such as ice crystals. Thus, an interferogram detection algorithm based on the vector pair method (GPV) was implemented. This method is very powerful both in terms of calculation times (13 times faster than CHT) and accuracy. It enables us to analyse each interferogram to deduce the LWC and IWC in the volume of measurement.

The frequency of the interferograms scattered by the droplets is measured using a 2D Fast Fourier Transform and the size of the crystals is assessed by measuring the width of the autocorrelation peak.

This processing coded in Matlab is performed at an average rate of 0.65 s per image (for 10 interferogram per image) on an office laptop. This only corresponds to three times as much as is necessary to load this image with Matlab (around 0.2 s). This computation time can be reduced by at least a factor 10 with optimised C++ programming.

Regarding planes' certifications rules, it is planned to use synthetic images to optimise the final design of the instrument, especially the defocus ratio.
A further step is to test the algorithm, as well as the
experimental device (currently undergoing optimisation), in
a calibrated icing wind tunnel (in terms of LWC and IWC).
The objective is to carry out flights test on an airplane with a
comparison with the Cloud Particle Imager or the Small Ice
Detector 3.

Acknowledgments

The research leading to this paper has received fund-
ing from the HAIC project, from the European Union's
Seventh Framework Program for research, technological
development and demonstration under grant agreement n°
ACP2-GA-2012-314314.

References

Adrian R J and Yao C S 1985 Pulsed laser technique application
to liquid and gaseous flows and the scattering power of seed
materials *Appl. Opt.* **24** 44–8

Brunel M and Shen H 2013 Design of ILIDS configurations for droplet characterization *Particuology* **11** 148–9

Brunel M, Coëtmellec S, Gréhan G and Shen H 2014a
Interferometric out-of-focus imaging simulator for irregular
rough particles *J. Eur. Opt. Soc.: Rapid Publ.* **9** 14008

Brunel M, Shen H, Coëtmellec S, Gréhan, G and Delobel T
2014b Determination of the size of irregular particles using
interferometric out-of-focus imaging *Int. J. Opt.* **2014** 143904

Brunel M, Ouldarbi L, Shen H, Gréhan G and Coëtmellec S 2014c Limits of rough particle sizing with low out-of-focus interferometric imaging *Particulology* **22** 52–7

Brunel M, González Ruiz S, Jacquot J and van Beeck J P A P 2015 On the morphology of irregular rough particles from the analysis of speckle-like interferometric out-of-focus images *Opt. Commun.* **338** 193

Chuang P Y, Saw E W, Small J D, Shaw R A, Sipperley C M, Payne G A and Bachalo W D 2008 Airborne phase Doppler interferometry for cloud microphysical measurements *Aerosol. Sci. Technol.* **42** 685

Cooley J W and Tukey J W 1965 An algorithm for machine
computation of complex Fourier series *Math. Comput.*
**19** 297

De Araujo Coelho A, Brenguier J L and Perrin T 2005 Droplet
spectra measurements with the FSSP-100: part 1. Coincidence
effects *J. Atmos. Ocean. Technol.* **22** 1756

Dehaeck S and van Beeck J P A P 2007 Designing a maximum
precision interferometric particle imaging set-up *Exp. Fluids*
**42** 767–81

Daszykowski M, Walczak B and Massart D L 2001 Looking
for natural patterns in data. Part 1: density based approach
*Chemom. Intell. Lab. Syst.* **56** 83

Debye P 1908 Das elektromagnetische feld um einen zylinder und
die theorie des regenbogens *Phys. Ztg.* **9** 775

Ester M, Kriegel H, Sander J and Xu X 1996 A density-based
algorithm for discovering clusters in large spatial databases

with noise *Proc. 2nd Int. Conf. on Knowledge Discovery and Data Mining (Portland, Oregon, 1996)* p 226

Glover A R, Skippon S M and Boyle R D 1995 Interferometric laser imaging for droplet sizing: a method for droplet-size measurement in sparse spray systems *Appl. Opt.* **34** 8409

Goodman J W 2009 *Speckle Phenomena in Optics. Theory and Applications* (Englewood: Roberts and Company Publishers)

König G, Anders K and Frohn A 1986 A new light-scattering technique to measure the diameter of periodically generated moving droplets *J. Aerosol. Sci.* **17** 157

Kimme C, Ballard D and Skalansky J 1975 Finding circles by an
array of accumulators *Proc. ACM* **18** 120–2

Lemaitre P, Porcheron E and Nuboer A 2007 Development of rainbow refractometry and out of focus imaging to characterize heat and mass transfers in two phase flow and aerosol collection by droplets *Proc. of ICONE 15, 15th Int. Conf. on Nuclear Engineering (Nagoya, Japan, 22–26 April 2007)*

Lawson R P, O'Connor D, Zmarzly P, Weaver K, Baker B, Mo Q and Jonsson H 2006 The 2D-S (stereo) probe: design and preliminary tests of a new airborne, high-speed, high-resolution particle imaging probe *J. Atmos. Ocean. Technol.* **23** 1462–77

Maeda M, Kawaguchi T, and Hishida K 2000 Novel interferometric measurement of size and velocity distributions of spherical particles in fluid flows *Meas. Sci. Technol.* **11** 13–5

Mounaïm-Rousselle C and Pajot O 1999 Droplet sizing by Mie scattering interferometry in a spark ignition engine *Part. Part. Syst. Charact.* **16** 160

Porcheron E, Lemaitre P, van Beeck J P A J, Vetrano R, Brunel M, Grehan G and Giraud L 2015 Development of a spectrometer for airborne measurement of droplet sizes in cloud *J. Eur. Opt. Soc.-Rapid 10* 15030–6

Pruppacher H R and Klett J D 2012 *Microphysics of Clouds and Precipitation* (Dordrecht: Springer)

Quérel A, Lemaitre P, Brunel M, Porcheron E and Gréhan G 2010
Real-time global interferometric laser imaging for the droplet
sizing (ILIDS) algorithm for airborne research *Meas. Sci.
Technol.* **21** 015306

Rad A, Faez K and Qaragozlou N 2003 Fast circle detection using
gradient pair vectors *Proc. VIth Digital Image Computing:
Techniques and Applications (Sydney)* ed C Sun et al

Shen H, Coëtmellec S, Gréhan G and Brunel M 2012
Interferometric laser imaging for droplet sizing revisited:
elaboration of transfer matrix models for the description of
complete systems *Appl. Opt.* **51** 5357

Shen H 2014 Interferometric out of focus imaging and digital
holography for the characterization of droplets or bubbles:
theory, optical design, characterization of a flow, synchronized
experiments *PhD Thesis* University of Rouen

Ulanowski Z, Hirst E, Kaye P H and Greenaway R S 2012
Retrieving the size of particles with rough and complex
surfaces from 2D scattering patterns *J. Quant. Spectrosc.
Radiat. Transfer* **113** 2457–7

Ulanowski Z, Kaye P H, Hirst E, Greenaway R S, Cotton R J, Hesse E and Collier C T 2014 Incidence of rough and irregular atmospheric ice particles from Small Ice Detector 3 measurements *Atmos. Chem. Phys.* **14** 1649

Van de Hulst H C and Wang R T 1991 Glare points *Appl. Opt.* **30** 4755–8

Wiener N 1930 Generalized harmonic analysis *Acta Math.* **55** 117