1.

Introduction

Speckle analysis and speckle interferometry have been used to obtain deformations, topography, and roughness characteristics of a wide variety of samples. Also, when considered dynamically, speckle helps to characterize the temporal evolution of changes in the shape of observed objects.¹ Speckle analysis techniques use the variations in intensity, polarization, and statistics to obtain the desired parameters on the sample. Most of the analysis techniques are heuristic in origin, and their objective is mainly practical and adapted for a given environment or application. Artificial intelligence techniques have been incorporated, along with learning algorithms (supervised and nonsupervised), to segment and identify regions of interest within data.^2,3 However, there is not a unique strategy for analyzing dynamic speckle images applicable to a wide variety of situations.

In this contribution, we merge together the outcomes of two well-established methods: principal component analysis (PCA) and image descriptors. PCA is a multivariate technique that has been successfully applied to characterize several kinds of spatial–temporal phenomena,⁴ such as $1/f$ noise,^5,6 digital compression algorithms,⁷ and dynamical images. Its statistical foundation makes this approach robust and reliable, providing spatial and temporal information of interest for the given case. Statistical image descriptors have been proposed and used to summarize in a single number or as a map, the overall characteristics of the phenomena.^8–10 These descriptors are usually simpler than the results from PCA and are computationally more economic from the point of view of time of calculation and resources allocation. On the other hand, PCA requires an additional layer of knowledge to be included within the analysis. As any other computational method, PCA always produces a result. This result needs expert interpretation and makes PCA more user-dependent, especially when reaching meaningful conclusions for a given problem. Only when incorporating the measurement conditions and the character of the data set, it is possible to link PCA results to meaningful physical magnitudes.

The application of PCA to biospeckle images has been already considered with quite promising results.^9,11,12 In this paper, we apply PCA to a given sequence of dynamic speckle images. PCA allows filtering sequences with specific spatial–temporal patterns. Then typical image descriptors are evaluated and compared with the classical dynamic speckle analysis methods applied to original and PCA-filtered sequences. Using this strategy, it is possible to determine which PCA band reveals the region of interest for each descriptor. Also, the time scales of the different processes can be determined. The input data used in this paper consist of three video sequences of dynamic speckle. One of the sequences corresponds with the drying evolution of a painted coin, and the other two are obtained by registering the speckle images of a bruised apple at different times.

Section 2 describes how PCA works for a temporally ordered collection of speckle frames (a dynamic speckle sequence). At the same time, we briefly present the results obtained from PCA and how they can explain the inner dynamics of the data. Then, in Sec. 3, we define the most useful image descriptors that are calculated for each case obtained before. Finally, Sec. 4 summarizes the main conclusions obtained in this contribution.

2.

Principal Component Analysis

PCA has been widely used in a variety of fields and situations. It has been applied to reduce the dimensions of data and to extract their most meaningful components and structures. Dynamic speckle sequences are collections of frames taken at a given temporal rate. The light source is coherent and speckle patterns change dynamically due to variations in the observed specimen. However, as far as the object is globally static along the measurement, frames show correlation among them. This correlation is analyzed and sliced by PCA. From a geometrical point of view, PCA can be seen as a rotation transformation from an $N$-dimensional coordinate system (the original frames) to a new $N$-dimensional coordinate system, where the new frames show no correlation. A dynamic speckle sequence has $N$ frames, $F_i$, each one having $M$ pixels organized as a rectangular image. Each rectangular frame is rearranged as a vector containing as many elements as pixels in one frame, $M$. This reshaping of data is fully reversible and rectangular frames can be retrieved. To properly evaluate the desired parameters, we have transformed each vector containing each frame to zero-mean vectors. This is easily done by subtracting the mean of each frame to each one, ${\overline{F}}_i=F_i-⟨F_i⟩$, where $⟨⟩$ denotes averaging. The $N$ vectors in the sequence are organized as an $N\times M$ matrix, and its covariance matrix, $S$, is calculated as an $N\times N$ matrix. The elements of the variance matrix of the original data set, $S_{ij}$, are calculated as the scalar product of the zero-mean vectors representing frames $i$ and $j$. Due to the character of the acquired dynamic speckle images, this matrix is not diagonal and shows correlations among original frames. PCA works by diagonalizing this covariance matrix, $S$.^4,5 The results obtained from PCA are $N$ eigenvalues (appearing at the diagonal of the new diagonalized matrix), ${\lambda }_{\alpha }$; $N$ principal components, ${\mathrm{PC}}_{\alpha }$ (uncorrleated frames); and $N$ eigenvectors, $E_{\alpha }$ (describing the transformation from original to-and-from principal component frames), where $\alpha$ runs from 1 to $N$ (we have chosen $\alpha$ to denote the index for the principal components and $i$ for the original frames). Eigenvalues describe the variance associated with each principal component. Most of the time it is more useful to evaluate the relative weight of each eigenvalue, $w_{\alpha }={\lambda }_{\alpha }/\sum _{\alpha =1}^N{\lambda }_{\alpha }$. The algorithm involved in the calculation of the principal components organizes eigenvalues in decreasing order (${\lambda }_{\alpha }>{\lambda }_{\alpha +1}$). Then the relative importance of a given principal component decreases when increasing its rank, $\alpha$. Principal components, ${\mathrm{PC}}_{\alpha }$, can be interpreted as speckle pseudoimages showing no correlation among them. Finally, eigenvectors can be organized as an $N\times N$ unitary matrix describing the rotation between the original frames, $F_i$, and the principal components, ${\mathrm{PC}}_{\alpha }$. In summary, PCA allows to move from a correlated set of frames (original speckle sequences) to an uncorrelated sequence (principal components).^4,5 This transformation is described as

Although principal components show no correlation, some subsets of them could be connected by inner data uncertainties. This fact allows grouping principal components into processes.⁵ The grouping strategy analyzes the uncertainties of the eigenvalues and links together those consecutive eigenvalues with overlapping uncertainties. The eigenvectors associated with isolated eigenvalues (not overlapping uncertainties) typically show a quasiharmonic variation associated to a temporal frequency.⁶ Considering this grouping method and its temporal behavior, it is possible to classify the set of $N$ principal components in classes showing different correlation properties and spatial–temporal patterns. This will be explained more in detail in Sec. 2.1. As a consequence of this classification strategy, it is possible to reconstruct a filtered version of the original sequence taking into account a selected collection of principal components. This is called principal component rectification, or PCA filtering, and it is given as a variation of Eq. (1)

In this paper, we have applied PCA to three collections of frames obtained from two types of samples. Figure 1 shows the average frame for each collection. This averaging is made along the temporal dimension, and it is different from the individual frame spatial averaging used to apply the PCA method. The samples were illuminated with an expanded and attenuated 10 mW HeNe laser. The speckle images were then registered by a CCD camera (Pulnix TM-6CN, CCIR Standard, interlaced scanning, pixel size $8.6\mu \mathrm{m}\times 8.3\mu \mathrm{m}$, objective lens $f=50\mathrm{mm}$, and $F^{\#}=16$), digitized by a frame grabber to 8 bits and stored in the memory of a personal computer.¹³ The first collection of data, collection A, is obtained from the drying process of a painted coin. In this case, a uniform film of paint was applied to an Argentinean 5 cents coin with topographic details in relief. The painted surface was then illuminated by the laser and 398 images ($512\times 512{\text{pixels}}^2$) were registered. From Fig. 1, the shape of the coin is clearly seen in collection A. Collections B and C are dynamic speckle sequences of an apple mechanically beaten. Collection B is obtained immediately after the beat and collection C is taken after half an hour. To detect the bruised region in the apple, we registered a sequence of 500 speckle images ($300\times 300{\text{pixels}}^2$) from a red delicious apple. The damage was caused by a controlled impact produced by letting fall a steel ball (diameter 21.9 mm and weight 133.6 g) from 20 cm height on the apple. The impact point is at the center of the bottom region of the image. For control purposes, a nonactive region (a steel blade) was included in the scene at the top-right corner of the image. The bruise could not be appreciated by visual inspection.

Fig. 1

Temporal mean of the frames for each data collection. Collection A is for the coin and collections B and C are for the apple.

2.1.

Results from the Principal Component Analysis

After applying the grouping and classification strategies described previously, we begin analyzing the eigenvalues and the temporal spectrum of the eigenvectors. In the first row of Fig. 2, we plotted the eigenvalues for each collection of frames in a semilog graph. The vertical lines limit three different classes classified using the PCA results.⁵ The first class in this classification, P1, comprises all the principal components, eigenvalues, and eigenvectors characterized by an isolated eigenvalue. The third class, P3, collects all the components that are represented by eigenvalues grouped together in one single process. In this last case, the uncertainties of consecutive eigenvalues overlap successively. Finally, the second class, P2, includes intermediate principal components between classes P1 and P3. A summary of the results of this classification is given in Table 1. Because the first principal component is the one contributing the most, we have quantified its relative weight, $w_1$, that is expressed as a percentage is 41.53%, 53.72%, and 56.16% of the total variance for the three collections of frames, A, B, and C, respectively. From Fig. 2, we can also check that the harmonic character of the eigenvectors is much clearer for the first eigenvectors (see middle row in Fig. 2) than for the eigenvectors corresponding to principal components of higher rank number.

Fig. 2

The columns correspond with the (a) coin and (b and c) apple sequences, respectively. The relative weights, $w_{\alpha }$, are plotted in the first row. The vertical lines denote the separation between different classes of dynamics classified in this paper. The first five eigenvectors are plotted in the central row. We may see how the dependence is quite close to a harmonic variation. Finally, the third row shows as a map the spectrum of the eigenvectors for the three data sets.

Table 1

Classification of the principal components into classes for the three data sets.

Classes	# range	Relative weight: ∑α∈Pwα
P1 (A, B, C)	[1:30], [1:36], [1:34]	0.8718, 0.7692, 0.7803
P2 (A, B, C)	[31:68], [37:127], [35:113]	0.0749, 0.0945, 0.0818
P3 (A, B, C)	[69:398], [128:500], [114:500]	0.0534, 0.1363, 0.1379

In addition to these preliminary results obtained from PCA, after classifying the components into classes and using Eq. (2), we may reconstruct the original data taking into account only those principal components included in each class (see # range in Table 1). In Fig. 3, we present the original object before painting (image a), and one of the frames for the original and the three filtered sequences (b, c, d, and e, respectively) for data collection A. These images show how poor the information that can be extracted from a single speckle image is. Only the profile of the coin is faintly perceived as a local diminished contrast. Then to obtain relevant information, every sequence can be treated statistically using image descriptors applied to the PCA outcomes. Some preliminary results are revealed when evaluating the mean, standard deviation (STD), and the signal-to-noise ratio (SNR) of the frames for the third class, P3, for the coin data set (see Fig. 4). SNR is defined as the quotient between mean and STD maps. From these figures, we can already see the hidden topography. STD is also considered as a descriptor in Sec. 3.

Fig. 3

The original object is represented when illuminated by (a) incoherent light and (b) coherent light. Images (c), (d), and (e) are frames extracted from the three filtered sequences P1, P2, and P3, respectively.

Fig. 4

The maps represent (a) the mean value and (b) STD for each pixel of the frames calculated along the whole sequence (398 frames) obtained for the third class (P3) of the coin data set. (c) SNR is calculated as the ratio between the mean and STD for each pixel. Hidden topography is revealed from the STD and signal-to-noise maps.

The main advantage of PCA is revealed when analyzing the temporal characteristics described by the corresponding spectra (see bottom row in Fig. 2). Actually, as we have seen before, the grouping strategy defines three classes of sequences: P1, P2, and P3. For class P1, there exists a linear relation between the frequency associated to each eigenvector and its rank.^6,14 Once this characteristic frequency is defined, it is possible to consider the relation between eigenvalues (denoting the importance of the principal component in the whole data set and given as a portion of the total variance of the data) and associated frequencies (obtained from the quasiharmonic dependence of the eigenvectors), as a power spectral density (PSD).⁶ Actually, when analyzing dynamic speckle sequences, this PSD follows a $1/f$ pattern, where $\mathrm{PSD}(\omega )=A_{1/f}/{\omega }^{\alpha }$, being $\alpha \simeq 1$.¹² For classes P2 and P3, it is not possible to define a single frequency for each eigenvector. To simplify the analysis, we extend to classes P2 and P3 the linear relation obtained for class P1 between the rank and the characteristic frequency. The frequency obtained from this linear trend is considered as a central characteristic frequency associated to each principal component (see Fig. 2). This frequency and its corresponding time constants are used to define the temporal and frequency limits of the classes identified previously. Table 2 shows the time scales for the three data sets. Class P1 describes dynamic behavior with frequencies between 0 and $f_{\lim ,1-2}$. Class P2 corresponds with frequencies between $f_{\lim ,1-2}$ and $f_{\lim ,2-3}$. Finally, class P3 expresses temporal frequencies above $f_{\lim ,2-3}$ until the Nyquist frequency. The same discussion can be given in terms of time scales. The values in Table 2 have been obtained knowing that the frame rate is 30 fps for the coin data set (collection A) and 1.25 fps for the apple data set (collections B and C).

Table 2

Limit frequency and characteristics time constant

Data set	flim,1−2 (Hz), τlim,1−2 (s)	flim,2−3 (Hz), τlim,2−3 (s)
A (coin)	1.125, 0.89	2.55, 0.39
B (apple $t=0\mathrm{h}$)	0.045, 22.2	0.159, 6.30
C (apple $t=0.5\mathrm{h}$)	0.043, 23.5	0.141, 7.08

From the previous analysis, it is possible to estimate the minimum resolvable spatial detail (spatial resolution) and the minimum resolvable time variation (temporal resolution). The spatial resolution is limited by the size of the speckle grain at the object plane, $G$, given as¹⁵

Eq. (3)

$$G\simeq 1.2\lambda F^{\#}(1+1/M),$$

where $M$ is the lateral magnification of the image-forming system. For the coin data set, $M=0.2$. Then, using Eq. (3) and the experimental conditions, the size of the speckle grain is $G\simeq 75\mu \mathrm{m}$. The estimation of the temporal resolution should consider the interaction between the $1/f$ PSD previously described¹² and other sources of noise.^6,16 A conservative estimation of the temporal evolution can be taken as limited by the characteristic time constant corresponding to the last frequency included within the first class. For the coin data set, this temporal limit would be ${\tau }_{\lim ,1-2}=0.89$ s (see Table 2). By combining the spatial and temporal resolutions, it is possible to define a speckle speed limit corresponding to the first class as $v_{\mathrm{limit}}=G/{\tau }_{\lim ,1-2}$.

3.

Image’s Descriptors

Global and local descriptors have been proposed to determine spatial and temporal characteristics of a dynamic speckle sequence. Some of the findings of the analysis shown in this paper are related to the spatial identification of temporal phenomena. Then it seems reasonable to focus our attention on time-domain descriptors. They are given in terms of the value of the pixel located at a position $(i,j)$, at a given frame, $k$, of the full $N$ frames sequence: $x_k(i,j)$. Because of its importance in medical applications, we also include laser contrast analysis (LASCA) descriptor that requires only one frame for its calculation. These descriptors will be easily represented as maps. They are listed as follows:

3.1.

Analysis of Image’s Descriptors Using Principal Component Analysis

Figure 5 shows the descriptors, previously defined, and calculated for the coin experiment (arranged in rows). The evaluation has been done for the original sequence of frames and the filtered sequences corresponding to P1, P2, and P3 classes (organized in columns). The same pseudocoloring palette was used in all of them to decrease subjective bias. This analysis demonstrates that only WGD, AVD, and Fujii’s difference were able to distinguish hidden topography from the original dynamic speckle coin sequence. However, when any of the previously considered image’s descriptors were evaluated for the PCA-filtered sequence using only class P3 principal components [$P=\{69,70,\dots ,398\}$ in Eq. (2)], all descriptors provide meaningful information of the hidden topography under the paint. In addition, none of these global parameters provided topographic information when using filtered sequences corresponding to class P1 [$P=\{1,2,\dots ,30\}$ in Eq. (2)]. Using the temporal parameters given in Table 2, we could conclude that dynamics above $f_{\lim ,2-3}=2.55\mathrm{Hz}$ reveals topography using any image descriptors. This means that high-frequency dynamics is dependent on the paint thicknesses, i.e., thinner regions are less active at high frequency than thicker areas. Contrarily, sequences involving frequencies below $f_{\lim ,1-2}=1.125\mathrm{Hz}$ do not give topographic information independent of the image’s descriptor used in the analysis. In this case, PCA has provided a quite simple way to decouple processes acting at different temporal scales. Even more, some of the results have provided important information about hidden features of the object under analysis, even when observed within a temporal window much smaller than the typical time scales of the physical mechanism under study.

Fig. 5

Maps for the descriptor applied to the original and filtered sequences for the coin data.

The bruised apple sequences have been treated in the same way as the coin. However, in this case, we are more interested in retrieving some information about the bioactivity of the sample. This is why we have performed the analysis with two sets of frames taken at different times: one sequence (collection B) is taken immediately after the apple has been bruised, $t=0\mathrm{h}$ and the second sequence is taken after half an hour, $t=0.5\mathrm{h}$ (collection C). In the analysis of the apple, we have focused our attention on the STD and Fujii descriptors because they showed a better performance when distinguishing the effect of the bruise. After evaluating them for original and filtered sequences of the two sets of data, they can be compared by calculating the difference in the descriptor’s map obtained at two different times after the beat. This analysis is shown in Fig. 6. The left column shows $\mathrm{Fujii}(t=0.5\mathrm{h})-\mathrm{Fujii}(t=0.0\mathrm{h})$, and the right column represents $\mathrm{STD}(t=0.5\mathrm{h})-\mathrm{STD}(t=0.0\mathrm{h})$. The location of the damage (center bottom of the image) is identified by this difference even when calculated with the original frame set (first row in Fig. 6). However, the capability of PCA to classify elements according to its temporal behavior allows a more detailed analysis. The beaten region is identified from filtered images P1, P2, and P3 (second, third, and fourth row in Fig. 6). Actually, although visible, this identification is a little harder for Fujii parameter calculated for the P1 class. With this analysis, we can infer some information about the dynamics of the process. Both Fujii and STD are larger at the beaten region just immediately after the beat than 30 min later for P2 and P3 classes. However, Fujii and STD are smaller at the beaten region for $t=0\mathrm{h}$ than for $t=0.5\mathrm{h}$. As these classes are linked with temporal frequencies, this differential behavior means that lower frequency dynamics are more important half an hour after the beat than immediately after the damage is produced. The opposite can be said for higher frequencies (larger than $f_{\lim ,1-2}$) denoting that the beaten region stops changing at high frequencies. In this discussion, it is important to note that the importance of each class (evaluated by the relative weights presented in Table 1) does not change significantly from collection B (just after beating the apple) and collection C (half an hour after the beat). This approximate constancy emphasizes the reliability of PCA to analyze spatial–temporal patterns.

Fig. 6

Maps of the difference between $t=0.5\mathrm{h}$ and $t=0\mathrm{h}$ for (a) Fujii and (b) STD image descriptors. This difference is evaluated for the original set of frames (first row) and for those obtained after grouping together the principal components into classes P1 (second row), P2 (third row), and P3 (fourth row). The temporal frequencies involved in these classes are higher for class P3 and are lower for class P1. Class P2 corresponds with intermediate frequency variations. The color bar applies to both parameters only by adding a $\times {10}^{-3}$ factor to the Fujii images.

4.

Conclusions

The application of PCA to dynamic speckle sequences obtained from industrial processes and biological samples makes possible a spatial–temporal description of the data. In this contribution, we demonstrate this capability by studying two cases that appear in the literature: the drying process of paint over a given topography and the bruising of a fruit. After applying PCA, we see that the most contributing images are associated with harmonic temporal variations. In addition, these principal components are independent and cannot be connected through the uncertainty of data. On the other hand, those principal components with the lowest contributions can be grouped together because of the inner statistics relations among data. This grouping strategy defines three classes of sequences. The first one represents most of the variance of the original data. The third class represents noise variations and rapid fluctuations. In between these two classes, there is another one that defines intermediate time scales. Temporal frequencies are associated with each class. This association is clearer for the most contributing class corresponding with low-frequency variations. Even more, in this case, it is possible to define a power spectrum density.

For the coin drying data, the analysis of the original data reveals the hidden topography for the WGD, AVD, and Fujii parameters. However, this analysis cannot be associated with the quantitative evaluation of the involved time scales until studying the PCA results, where the three filtered sequences are associated to three different time scales and frequencies. In the coin case, the application of the descriptors shows that the topography under paint is mostly contained in class P3, corresponding to low eigenvalues. For the bruised apple experiment, the result shows that both the STD and Fujii descriptors identify the damaged region even from the original sequences. Again, only after analyzing PCA results and comparing them using sequences taken at different times after the bruising, is it possible to perceive some changes. We have seen that low-frequency dynamics is more important 30 min after the beat at the location of the damage. The opposite happens with medium- and high-frequency contributions. This means that half an hour after the beating, those processes involving rapid variations decrease.

Summarizing these results, we conclude that PCA is capable of revealing hidden contributions to the speckle variations. These contributions are easily related to separable time scales. This fact makes possible a quantitative and qualitative spatial–temporal description of the phenomena under analysis. In addition, the application of speckle descriptors to filtered sets of data improves the capability of identification of underlying biological or industrial time-dependent processes. Further improvement can be expected, if different descriptors are combined by using supervised or no supervised neural networks.