Dynamic speckle analysis using multivariate techniques

Dynamic speckle patterns inform about the temporal characteristics of the underlying physical phenomena in the sample. Usually, several phenomena concur and contribute to the final observed intensity variations, including scattering, Doppler effect, polarization rotation, etc [1]. Since the same effects can be due to one or more of these origins, it is difficult to identify the actual mechanisms at work, and the quantification of their individual contributions.

Several techniques for analyzing dynamic speckle patterns have been proposed in the literature [1]. Many of them are heuristic in origin and others are specific to some applications. Recently, some artificial intelligence tools have been developed. They include both supervised [2] and non-supervised learning algorithms [3]. They make use, simultaneously, of several measuring techniques for screening some effects. Nevertheless, the problem is not fully solved and the identification and segmentation of regions showing different inner dynamics is of considerable interest, for example, in biology and agronomy. In this work we analyze sequences of images of dynamic speckle using a multivariate technique: the principal components analysis (PCA) [4]. This method can distinguish processes occurring in a video sequence which is characterized by different correlation structures. The computational resources and time necessary to perform PCA are proportional to the data volume (number of frames × number of pixels at each frame) and it is in the range of tens of seconds. This processing time is much longer than the time consumed in evaluating other figures of merit [5]. However, PCA retrieves more information about the dynamic characteristics of the physical process providing temporal parameters beyond the acquisition and processing time duration, along with spatial structures. To show how PCA works in dynamic speckle analysis we apply it to a quite explored experimental situation: the drying of paint on a surface of hidden topography [1, 6]. In this case the sample is a coin. This kind of object has been used in the literature to compare the performance of different algorithms.

To understand the method and its limitations, in section 2 we briefly present the PCA technique and its practical implementation for the analysis of a video sequence. Section 3 explains how the results of the PCA describe the inner dynamic of the physical process that is observed, and how the results can be applied to any other video sequence of dynamic speckle patterns. Finally, in section 4 we summarize the main findings of our contribution.

PCA, a multivariate technique, analyzes the structure of the variance of a collection of statistical realizations of a set of variables [4]. Among other applications, it reduces the dimension of a set of data by generating a collection of new variables that do not exhibit correlation among them, and selecting the most meaningful. These new variables are linear combinations of the original, and therefore they represent the same magnitude. In this work, the original variables are the collection of N consecutive frames organized as a video sequence, grabbed by a CCD camera from an active sample producing dynamic speckle pattern when illuminated by laser light. In this paper, we name F the following N-dimensional variable: $F=\left\{ {{F}_{1}},{{F}_{2}},\ldots ,{{F}_{t}},\ldots ,{{F}_{N}} \right\}$, where F_t represents the t frame in the temporal sequence. To properly account for inner correlations between frames, the mean of each frame is removed from each one. The PCA transformation provides a new set of frames that we label as $\bar{F}$. We also know that each frame represents a bi-dimensional distribution of irradiance having R rows and C columns. From an statistical point of view, each pixel is a realization of F. Then, the number of realizations of the original variable F is $M=R\times C$. If we maintain this bi-dimensional arrangement, we would need to deal with a data cube of consecutive frames. However, it is possible to transform each R × C frame into a M vector by a simple rearrangement of the data. Both the removal of the mean and the rearrangement of the pixel's values are reversible procedures, so the original frames can always be retrieved.

The frames are connected through the co-variance matrix, S, of the zero-mean original data set, $\bar{F}$. As expected, S shows non diagonal elements accounting for the inner correlation among frames. However, it is possible to diagonalize this co-variance matrix by moving to another set of variables that now are uncorrelated. This transformation, that can be interpreted as a rotation within a N-dimensional space, is made by the PCA method. This rotation moves from the frame reference system, F, where the images are correlated, to a new reference system, PC, where the correlation between frames is zero. Figure 1 shows the transformation, including the original frames (${{F}_{1}},{{F}_{2}}$, and F₃) that can be interpreted along directions $(1,0,0,\ldots ,0)$, $(0,1,0,\ldots ,0)$, and $(0,0,1,\ldots ,0)$, respectively. The new uncorrelated variables, named as principal components, are related to the original variables by the following linear relation:

$$\begin{eqnarray}&&{\rm P}{{{\rm C}}_{\alpha }}=\mathop{\sum }\limits_{t=1}^{N}{{e}_{\alpha ,t}}{{\bar{F}}_{t}}.\end{eqnarray} \tag{ 1 }$$

The dimensionless coefficients, ${{e}_{\alpha ,t}}$, of the linear combination describing this transformation are the eigenvectors, ${{E}_{\alpha }}$, of S. In figure 1 we interpret these coefficients as the new directions, e₁, e₂ and e₃, along the new PC reference system. In this figure we have also presented the images corresponding to the given principal components. Although the original images—represented as F₁, F₂, F₃—are quite similar among them, the new images—described by the principal components PC₁, PC₂, and PC₃—are already different and show no correlation. The criterion to obtain each direction is that the corresponding principal component should contain the maximum amount of variance of the remaining data set. As a matter of fact, the variance of the resulting component images will be an important result of the analysis, and the method provides these values in decreasing order. It can be demonstrated that the ${{e}_{\alpha ,t}}$ coefficients obey the following relation:

$$\begin{eqnarray}&&\left( S-{{\lambda }_{\alpha }}I \right){{E}_{\alpha }}=0,\end{eqnarray} \tag{ 2 }$$

where I is the $N\times N$ identity matrix. The eigenvalues ${{\lambda }_{\alpha }}$ obtained from the previous equation represent the variance associated with each principal component, ${\rm P}{{{\rm C}}_{\alpha }}$. Therefore, the fraction of the total variance explained by each PC, ${{\Omega }_{\alpha }}$, is given by

$$\begin{eqnarray}&&{{\Omega }_{\alpha }}=\frac{{{\lambda }_{\alpha }}}{\mathop{\sum }\limits_{\alpha =1}^{N}{{\lambda }_{\alpha }}}.\end{eqnarray} \tag{ 3 }$$

We may also note that the principal components are organized in decreasing order of their contribution to the total variance of the data, being the sub-index $\alpha =1$ the one corresponding with the most relevant PC. The eigenvectors, ${{E}_{\alpha }}$, can be arranged as a $N\times N$ matrix. If both the original variables, ${{\bar{F}}_{t}}$, and the principal components, ${\rm P}{{{\rm C}}_{\alpha }}$ are organized as two $M\times N$ matrices, the transformation given in equation (1) can be written as: ${\rm PC}=\bar{F}E$. Using this matrix relation, or equation (1), we obtain the original frames from the principal components as:

$$\begin{eqnarray}&&{{\bar{F}}_{t}}=\mathop{\sum }\limits_{\alpha =1}^{N}{{e}_{\alpha ,t}}{\rm P}{{{\rm C}}_{\alpha }}.\end{eqnarray} \tag{ 4 }$$

This equation reconstruct the original data (except for the mean of each frame that can be easily added to the result of equation (4)) using all the principal components. However, the sum can be carried out by selecting some meaningful principal components to produce a filtered, or rectified, version of the original data.

Zoom In Zoom Out Reset image size

As a general tool of multivariate analysis, one of the main drawbacks of the PCA method is the accurate interpretation of the results for a given application. However, this statistical origin is also an advantage of the method because its sound foundation. In fact, PCA has been successfully applied to a variety of fields both in optics: spectroscopy, facial recognition, etc; and other areas as econometrics, or computational linguistics [4]. Our case is similar to the analysis of the noise in a sequence of video frames.

2.1. Grouping and classification of principal components

In this contribution PCA is applied to a sequence of 398 images obtained at a frame rate of 30 fps (sampling frequency, ${{f}_{{\rm Sampling}}}=30$ Hz) with a size of 512 × 512 = 262144 pixels and lasting 13.3 s. The analysis has been made with a personal computer with an Intel-i7 processor, 16 Gb of RAM, under the Matlab environment. The PCA results are obtained after 21 s. If the data follow a normal distribution, the grouping and classification are straightforward and consumes less than 1 s. However, if the normality assumption fails, then it is necessary to calculate fourth order cumulants (see reference [8]) and the grouping process takes around 4 min. Finally, the reconstruction of the frames using a selected set of frames is evaluated in 2 s. The computational time can be reduced when using dedicated and compiled algorithms.

Our sequence of images contains the dynamic speckle produced by a coin covered with paint along the drying of the paint coating [1]. The underlying topography cannot be perceived during the drying process, but the coating thickness varies accordingly, and thin paint dries before thick ones. So, measures of dynamic activity should show the hidden topography. The method can be straightforwardly adapted to other dynamic processes.

We first represent the relative importance of the PC to the total variance of the original data set (see figure 2). The first principal component explain 41.5% of the total variance of the data, and the following ones are contributing in decreasing order.

Zoom In Zoom Out Reset image size

Figure 3 shows a log–log plot of the percentage of the eigenvalue contributions obtained after diagonalizing the covariance matrix, S. They are represented as a function of the ordinal of the associated principal component. After analyzing the results, three classes of structures can be observed that we label as P1, P2, and P3.

Zoom In Zoom Out Reset image size

The first one, P1, is formed by the first 30 PCs. This class represents 87.2% of the total variance of the original data. These PCs are isolated and independent because they are single-component processes. This means that their eigenvalues, considering their uncertainty due to the sample, do not overlap. Their associated eigenvectors can be fitted to sinusoidal functions of time (see inset of eigenvector #6 in figure 4). Therefore, for each principal component included in this spatial-temporal structure, the associated variance , i.e., the eigenvalue, is related to a single frequency. This frequency is determined by the time evolution of the associated eigenvector [10, 11]. To prove this, figure 4 shows the Fourier transform of each eigenvector. We can observe that the sinusoidal eigenvectors included in the P1 class present a narrow peak if compared with the rest that are not isolated eigenvectors. The results obtained from the PCA for this first structure P1, can be seen as a PSD function [10, 11]. This interpretation is graphically presented in figure 5, that is a log–log plot of the eigenvalues. Versus the frequencies associated with each eigenvector. These data can be linearly fitted showing a temporal evolution of the type ${{f}^{-\alpha }}$. Then, the ${\rm PSD}(f)$ is described by the following equation:

$$\begin{eqnarray}&&{\rm PSD}(f)=\frac{30\pm 2}{{{f}^{1.09\pm 0.02}}}.\end{eqnarray} \tag{ 5 }$$

This fitting yields a value of $\alpha =1.09\pm 0.02$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The third class, P3, groups together in the same multi-component process principal components ranging from #70 to #398. They represent 5.4% of the total variance of the original data. In figure 4 we can see how the range in temporal frequency components of these PCs is wide and contains high frequency components until the Nyquist frequency (${{f}_{{\rm Nyquist}}}=0.5{{f}_{{\rm Sampling}}}=15$Hz).

Finally, the second class, P2, includes those PCs from #31 to #69 that cannot be included in P1 (representing a PSD) nor the ones in the P3 class (representing a uncorrelated high-frequency temporal noise). By analyzing the uncertainties in the data represented in figure 3 and applying the grouping strategy, we conclude that this subset of principal components is not a single process. Nevertheless, the eigenvectors of the P2 class share the feature of being defined by two frequencies or even three (see figure 4 and the inset representing eigenvalue #50), instead of one as it happens in P1 class.

Once the principal components are split according to the aforementioned criteria, we can proceed to the rectification, or filtering, process. Rectification consists in the reconstruction of a sequence of N frames (the same number as in the original sequence) using only the principal components obtained in each one of the defined classes. We select a set of principal components when performing the sum in equation (4) and add the mean value of each frame. We obtain three sequences of frames, one for each class. The importance of each class can be given by the sum of the contributions, ${{\Omega }_{\alpha }}$, of the principal components belonging to each one. The first class, P1, involves 87.2% of the total variance of the original data, P2 represents about 7.5%, and the third class, P3, accounts for remaining 5.4%.

Several types of analysis can be performed on the images to find the reason behind their differences (see for example Fujii's method in [1]). In figure 6 we show the mean and the standard deviation of each rectified frame set. For the three sequences, the mean of each pixel, and its standard deviation are calculated and shown as images. We see that the images of the mean do not reveal useful information. However, the standard deviation images reveal the underlying topography of the object. In this way, a high value in the image implies a point where the dynamic speckle varies with high amplitude and fast temporal evolution. A quite similar analysis, using the noise instead of the signal as the figure of merit, has been successfully applied to the detection and identification of buried land-mine objects with principal components [16].

Zoom In Zoom Out Reset image size

The first class, P1, shows a high absolute value of the standard deviation with an almost uniform spatial distribution over the image. In the second one, P2, regional differences appear, which are clearly visible in the third class, P3, where the topography of the coin, letters and numbers, can be distinguished. The region surrounding the digit '5' has the strongest activity. The fundamental difference between the second and the third classes is spatially located around the digit, clearly visible in the third but not so much in the second. The standard deviation decreases from P1 class values to P2 and P3. However, in these last two classes, the spatial distribution of the variance changes more than in P1, revealing the underlying topography. This analysis of the noise is possible after applying the PCA method and the previously explained grouping strategy.

PCA can also assign time scales to these classes. >From figures 3 and 4 we infer that the highest frequency achieved by class P1 corresponds to about 1 s (1 Hz), and about 0.3 s for class P2. This implies that the drying dynamics that produces P1 involves movements in the speckle in the order of seconds or longer, and P2 between 1 and 0.3 s. The third class, P3, is then, associated to an even shorter time scale.

We can relate topography and time sales with the drying mechanism. The drying of paint is a rather complex process involving several steps including solvent evaporation and refractive index variations, features to which dynamic speckle is sensitive [17]. The initial stage of the process is when the solvent evaporates, and it is called the constant rate period. Afterwards, there is a falling rate period in which the drying rate decreases to zero. The beginning of the falling rate period depends on the film thickness and solvent evaporation [18]. Since the thickness of the coating on the over-relief of the coin is locally variable, the paint can be locally in different drying stages. Thin paint regions should be in later stages compared to thick ones. Therefore, activity in the speckle pattern varies with the hidden coin topography [6]. The results obtained from the PCA method shows that classes P2 and P3, associated with time periods shorter than 1 s, reveal the hidden topography of the coin and can be related with regions at later drying stages. From the results obtained for the P3 class, we observer that the standard deviation is smaller in those regions where the paint coating in thinner. This means that the contribution to the rapid variation of the P3 class is proportional to the thickness of the paint. At the same time, class P1, is identified as having a time period longer than 1 s and following a $1/f$ PSD. We should emphasize that this information has been obtained from a sequence frame lasting 13.3 s. Although this time period is much shorter than the total drying time, the results obtained from PCA already reveal different processes that last longer than acquisition time.

We presented the application of a multivariate technique, the PCA method, to the analysis of a sequence of images of dynamic speckle interferometry. We obtained quantitative results of different processes that can be related with the inner dynamic of the drying process of paint coating a coin. The classes identified here are related with temporal and spatial distributions at different spatial-temporal scales. Additionally, these results have been obtained from a sequence lasting just 13.3 s. This time duration is much shorter than the typical drying times of paint, allowing fast monitoring of the process. The method provides a quantitative evaluation of the reported classes of spatial-temporal evolutions. Most of the variance of the original data (87.2%) is explained as a $1/f$ PSD with a time scale longer than 1 s. The remaining variance (13% of the total) presents a temporal evolution with time constants faster than 1 s. Furthermore, they already reveal the hidden topography that is associated with later drying states. The capability of the PCA method to section the variance of the data into independent components allows a meaningful analysis relating physical sources to different classes of spatial-temporal evolution.

The computational time duration of the PCA method makes difficult its real-time application with current and standard computational capabilities. In these conditions, PCA can be seen as a fast pre- and post-processing tool that allows filtering of the data and provides parameters with time characteristics longer than the time duration of the analyzed sequence. Also, PCA results are given in terms of PSDs and spatial maps of characterizing parameters. This identification relates the observed classes with the dynamic processes associated to evaporation of solvent, displacement of pigments at different sizes or scales. The method provides new information, allows the trimming of the components of interest, and does not require previous knowledge of the dynamics of the phenomenon under investigation as it adapts itself to the signal.