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

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.


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 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 (2014aBrunel et al ( , 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 −3 ). 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.

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 andvan 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 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)). (1)

( ) ( )
In this equation, θ is the scattering angle (figure 1), m is the refraction index, and λ 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 (φ) 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 . 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 (2014aBrunel et al ( , 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 (2014aBrunel et al ( , 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 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.

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 (θ) of 90°. The sensor used for this assembly is an sCMOS sensor (2560 × 2160 pixels 2 , 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 (α) 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 + Δp (figure 3). M L is the transfer matrix modelling the refraction by the lens. The total transfer matrix ( ) M 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)).
(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.

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 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.

Detection of each interferogram in the image
Detecting each interferogram (interference disk) in an outof-focus image consists in identifying each disk in the image   (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 ( → V j and → 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: 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: Therefore, for each gradient vector → V j of the image, a vector → 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 A k 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.

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 nonoverlapping areas of each interferogram, and the bottom line presents the 2D Fourier transforms of those interferograms. 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.

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.

Analysis of the interferograms scattered by drop-
lets. 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.

Analysis of the interferogram scattered by ice crys-
tals. For interferograms scattered by irregular particles like ice crystals, Brunel et al (2014b) suggested a simple relationship linking the speckle grain size ( ) δ to the dimension of a particle ( ) ∆ .
λ Β δ ∆= tot (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 ( ) δ 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 f is the power spectral density and 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 δ 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)). Likewise, it is simple to measure the orientation of the ice crystal by measuring the angle between its major axis and  the horizontal axis. To assess the IWC, the average speckle grain size ( ) δ 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.
In equation (8) S grain is the surface of the speckle grain. Then, the average diameter of the crystal ( ) ∆ is calculated using the equation below.
Thus, knowing the dimensions of the volume of measurement of each image ( = × × V 10 cm 10 cm 1 mm) and the size of the droplets (d ) and ice crystals (Δ), the IWC and the LWC are determined (equation (10) In equation (10) N image is the number of images on which the measurement is performed. 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.

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).
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(Brunel et al , 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 (Δ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.  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 B tot coefficient of the total transfer matrix is needed. B tot is deduced from equation (2): And z2 is computed with the conjugation relation: 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 LWC algo and IWC algo are calculated (equation (13)) and then compared to LWC syn im and IWC syn im computed from the input data of same synthetic images.
In these equations d i and ∆ 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 1000 image ), 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 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-offocus imaging) to measure the size of ice crystals.

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 (2014aBrunel et al ( , 2014bBrunel et al ( , 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.