Performance of intensity-based non-normalized pointwise algorithms in dynamic speckle analysis

E. Stoykova; D. Nazarova; N. Berberova; A. Gotchev

doi:10.1364/OE.23.025128

1. Introduction

Optical metrology with dynamic laser speckle analysis is gaining popularity due to widening scope of applications and simplicity of data acquisition [1]. This technique requires coherent illumination of a diffusely reflecting object and capture of speckle patterns formed on its surface. Any microscopic changes of the surface over time reflect into changes of the speckle patterns [2]. This makes possible non-destructive detection of physical or biological activity through statistical processing of the acquired speckle data [3,4]. Dynamic laser speckle has been applied to study blood flow perfusion in human tissues [5,6], bacterial response [7], plant growing processes [8,9], seeds viability [10,11], animal reproduction [12], drying of paints and coatings [13], fruits quality [14,15] and bread cooling [16]. Usage of modern two-dimensional (2D) optical sensors and powerful computers promoted development of algorithms, which rely on acquisition of a temporal sequence of correlated speckle images. Among them, the most perspective are the algorithms whose output is a 2D spatial contour map of the estimate of a given statistical measure. The estimate at each point of the map is built by using a time sequence of intensity values taken at one and the same pixel in the acquired speckle images. The main benefit from pointwise processing is visualization of varying activity across the object surface.

Speckle nature of the acquired patterns exhibits itself as a strongly fluctuating intensity from point to point. At a finite number of patterns, this leads to substantial fluctuations in space of any pointwise intensity-based estimate even for spatially invariable activity. These fluctuations inevitably worsen sensitivity and spatial resolution of the dynamic speckle method because they decrease the contrast of the built 2D activity map. In general, the spread of fluctuations at a given activity is affected by speckle statistical properties and should decrease with the length of the processed time sequence of intensity values that is used to build the estimate at a point. For processing one and the same sequence of speckle images, the spread crucially depends on the chosen algorithm. Therefore, the algorithm performance should be assessed through analysis of the statistical measure it calculates at a given activity. Another important parameter when comparing the pointwise algorithms is the computation time they need to build the activity map.

Up to now, quite a few papers discussed and compared performance and efficiency of different algorithms in dynamic speckle analysis. Analysis of reliability of a time history of speckle pattern for activity characterization is done in [17]. The time history is an image formed by taking an identical single column from each image in a temporal sequence of images and thus does not yield an activity map for the whole object surface. A protocol for a quality test of the acquired speckle patterns is developed in [18]. A comprehensive review of using frequency and space-frequency representations to describe activity is given in [19]. Comparison of quality of 2D activity maps produced by pointwise intensity-based algorithms composed as modifications of the generalized difference (GD) approach is made in [20] by applying them to process data from a paint drying experiment. In [21] analysis of 2D activity maps is done for pointwise algorithms based on evaluation of temporal correlation and structure functions but again by qualitative comparison of the maps obtained by processing of both synthetic and experimental data.

In the present paper we propose to base the performance assessment of pointwise intensity-based algorithms on the probability density function (PDF) of the estimates they build. We show, as a first task, that this statistical approach allows for quantitative evaluation of the algorithm sensitivity to activity variation across the tested object. As a second task, we evaluate performance of the pointwise intensity-based non-normalized algorithms which rely on the differences between intensity values at two moments. We compare activity estimation by the GD algorithm [22] and its improved modification with a running temporal window to the recently proposed by us usage of a structure function [21]. As a third task of the paper, we propose a new algorithm and prove that it provides a high quality activity map at less computation time. The analysis and comparison of the studied algorithms are done by computer simulation and their performance is illustrated by processing experimental data.

2. Pointwise intensity-based algorithms without normalization

In the experimental set-up for acquisition of speckle patterns in Fig. 1 (a) an expanded laser beam illuminates the object to be tested. A CCD or CMOS camera is adjusted to focus the object. The optical axis of the camera is normal to the object surface. The set-up is positioned on a vibration-insulated table. The camera with a pixel interval $Δ$ records successively a sequence of $N$ correlated images of size $N_{x} \times N_{y}$ for time $T$ at a sampling rate $1 / Δ t = T / N$ . A time sequence of 8-bit encoded intensities $I_{k l, n} \equiv I (k Δ, l Δ, n Δ t), n = 1.. N$ is formed at each point $(k Δ, l Δ), k = 1.. N_{x}, l = 1.. N_{y}$ of the acquired images as is shown in Fig. 1 (b). This data can be used to build a pointwise estimate of a given statistical measure over $T$ . The 2D spatial distribution of the estimate gives a map of the object activity.

Fig. 1 Experimental set-up for acquisition of a sequence of N speckle patterns for time $T$ for dynamic speckle analysis (a); forming a time sequence of intensity values at each point of the acquired patterns (b).

Download Full Size | PDF

We analyze the case of uniform intensity distribution within the illuminating laser beam and equal reflectivity all over the object. This means that the mean value and the variance, $σ^{2}$ , of intensity fluctuations remain the same across the object. We assume that activity observed on the object surface is developing as a result of a stationary process. Then the intensity fluctuations, $I_{k l, n}$ , are characterized at each pixel $(k Δ, l Δ), k = 1.. N_{x}, l = 1.. N_{y}$ by a temporal autocorrelation function (ACF), $R_{k l} (τ = m Δ t)$ , where $τ = m Δ t$ is the time lag and $m \geq 0$ takes integer values. The undergoing activity may have different time scales across the object surface. Statistically this variation can be described by the 2D spatial distribution $τ_{c} = τ_{c} (k Δ, l Δ)$ of the temporal correlation radius of $R_{k l} (τ = m Δ t)$ . At $τ > τ_{c}$ correlation between the intensity values is considered as weak. Hence, if the ACF at a pixel $(k Δ, l Δ)$ is given as a product $R_{k l} (τ) = σ^{2} ρ_{k l} (τ)$ of the variance, $σ^{2}$ , and the normalized ACF, $ρ_{k l} (τ)$ , at this point, the value of $τ_{c} (k Δ, l Δ)$ is defined as $ρ_{k l} (τ_{c}) = ρ_{0} \leq 0.5$ . The choice of $ρ_{0}$ is task-specific. The averaging time interval T is a crucial parameter for pointwise processing of correlated speckle patterns. Obviously, its length should be large enough to ensure a reliable statistical estimate and short enough to follow the change of activity in time

A popular characterization of activity is to build a GD-estimate [20,22] at a point $(k, l) \equiv (k Δ, l Δ), k = 1.. N_{x}, l = 1.. N_{y}$ as follows:

\hat{G} (k, l) = \sum_{i = 1}^{N} \sum_{j = i + 1}^{N} | I_{k l, i} - I_{k l, j} |,

Alternative forms of the estimate (1) and its relation to the estimate of the variance of intensity fluctuations at each point have been studied in [20,22]. As is shown in [22], this statistical measure accumulates differences formed between any two intensity values in the temporal sequence at

(k Δ, l Δ)

. For this reason the GD estimate is not able to provide a high-contrast activity map and exhibits low sensitivity. The algorithm performance is substantially improved if accumulation of differences is done within a window of length

t_{w} = w Δ t

as it has been proposed in [20]:

{\hat{G}}_{w} (k, l, w) = \sum_{i = 1}^{N - w} \sum_{j = i + 1}^{i + w} | I_{k l, i} - I_{k l, j} | .

Actually, some weight from 0 to 1 may be attached to each entry inside the window however we restrict our analysis only to the case when all weights are taken equal to unity. The estimate given by Eq. (2) was called a mean windowed difference in [20] but we prefer to use a modified generalized difference (MGD). The selectivity introduced by the window leads to much higher sensitivity of this algorithm.

Recently we introduced a temporal structure function (SF) as an estimator of activity [21]. For activity induced by a stationary process and a uniform spread of intensity fluctuations across the object, the SF at each spatial point is given by $S_{k l} (τ) = 2 σ^{2} [1 - ρ_{k l} (τ)]$ . The estimate of the temporal SF at a time lag $τ = m Δ t$ is determined by

\hat{S} (k, l, m) = \sum_{i = 1}^{N - m} {(I_{k l, i} - I_{k l, i + m})}^{2} .

We propose in the paper a new measure of activity by replacing the square in Eq. (3) with the absolute value of the difference

I_{k l, i} - I_{k l, i + m}

. The estimate is calculated as

{\hat{S}}_{m} (k, l, m) = \sum_{i = 1}^{N - m} | I_{k l, i} - I_{k l, i + m} |

and will be called a modified structure function (MSF) throughout the paper. The algorithms given by Eqs. (3)-(4) include only one summation and therefore can perform faster than the algorithms given by Eqs. (1)-(2). Actually the SF and MSF algorithms should be compared to

{\hat{G}}_{w} (k, l, w)

because, similarly to the window

t_{w} = w Δ t

, the time lag,

τ

, introduces selectivity and thus improves the contrast. They could replace the modified GD only if they guarantee the same quality of the activity map.

An activity increase leads to a higher value for each of the estimates (1)-(4). As the value itself it not meaningful (with the exception of the structure function which is related to the normalized correlation function), the thing that matters is the sensitivity of each algorithm, i.e. the smallest recognizable increment in activity, or the contrast of the provided activity map. As it can be seen, the estimates (1)-(4) are not normalized. They are composed from the differences of intensities observed at different moments at a given spatial point. Since the spread of intensity fluctuations in a speckle pattern depends on the laser beam intensity, normalization is required to cope with experimental data obtained at non-uniform illumination and varying reflectivity across the object surface. For example, the Fujii estimate [23] evaluates the mean value of the ratio between the absolute difference of intensities at neighboring pixels, $| I_{k l, i} - I_{k l, i + 1} |$ , and their sum $I_{k l, i} + I_{k l, i + 1}$ . The correlation-based estimates in [21] undergo pointwise normalization by using the estimate of the variance at each point [21]. The main advantage of the non-normalized estimates is that they provide activity maps with a higher contrast than the normalized estimates.

3. Histogram of the estimate fluctuations

Due to the limited length of T and the speckled nature of the recorded images, all pointwise estimates strongly fluctuate from point to point across the activity map. Fluctuations decrease its contrast and hence sensitivity and spatial resolution of each algorithm. Therefore it is crucial to characterize the algorithms through statistical description of the estimates they build.

We propose to evaluate performance of an algorithm by analyzing the PDF of the estimate it produces. For the pointwise algorithms, such evaluation can be done by processing a spatial region corresponding to one and the same activity in speckle images obtained for a synthetic or real object. The region should contain enough points to provide a reliable distribution of frequency counts or a histogram of the estimate values that can accurately represent its PDF. The task is alleviated by the low spatial correlation of the estimate values within the spatial region of the same activity. Spatial correlation is determined by the average speckle size on the object surface.

Usage of simulated data offers the advantage to construct an object with controllable variation of activity. So, we conducted a numerical experiment with a specially designed synthetic object composed by four regions of constant activity in order to establish the properties of the estimates depending on the ratio $T / τ_{c}$ . Actually, this ratio is a qualitative measure of information capacity provided by a given time sequence. The capture of speckle patterns for this object was simulated for a He-Ne laser emitting at 632.8 nm at uniform illumination and equal reflectivity. The simulation included the following steps:

1) generation of initial 2D array of random delta-correlated phase values uniformly distributed from 0 to 2 $π$ ;
2) generation of a 2D array $τ_{c} = τ_{c} (k δ, l δ)$ , $k = 1...2 N_{x}, l = 1...2 N_{y}$ with $δ = Δ / 2$ to produce varying activity across the sample;
3) generation of a sequence of 2D delta-correlated in space random phase distributions $ϕ (k δ, l δ, i Δ t), k = 1..2 N_{x}, l = 1..2 N_{y}, i = 1.. N$ , correlated in time by the normalized ACF $ρ_{k l} (τ) = \exp [- τ / τ_{c} (k, l)]$ where phase evolution was introduced as is described in Ref [24];
4) simulation of subjective speckle intensity at the optical sensor aperture; the complex amplitude of the light reflected from the object is $U_{S} = \sqrt{I_{0}} \exp {- j ϕ (k δ, l δ, i Δ t)}$ at the moment $i Δ t$ for the illuminating beam with intensity $I_{0}$ . In paraxial approximation the complex amplitude of the light falling on the CCD array is $U_{C C D} = F T^{- 1} {H \cdot F T (U_{S})},$
where H is the coherent transfer function of the registration system and FT[..] denotes the Fourier transform (for a diffraction limited 4f system H is reduced to a circ function [25]);
5) simulation of integration by the CCD camera pixels with size $Δ$ ; realizations ${| U_{C C D} |}^{2}$ were modeled as arrays of 512 × 512 pixels and their values were summed in a window of size 2 × 2 pixels to simulate the integration by the CCD camera.

The obtained arrays of $N_{x} \times N_{y}$ = 256 × 256 pixels were saved as bitmap images with 256 gray levels and processed by the studied algorithms. To construct the test object, the simulated dynamic speckle patterns were divided into four rectangular regions Z₁, Z₂, Z₃ and Z₄ of size 64 × 256 pixels each with $τ_{c}$ taking values of 10 $Δ t$ , 20 $Δ t$ , 40 $Δ t$ and 80 $Δ t$ , i.e. the correlation radius was increasing in geometric progression. We processed image sequences generated for three speckle sizes at acquisition time $T_{1} = 64 Δ t$ , $T_{2} = 128 Δ t$ , $T_{3} = 256 Δ t$ and $T_{4} = 512 Δ t$ . For the sake of brevity, we describe the three simulated speckle types as speckle “A”, “B” and “C” . The speckle size was derived from the 2D distribution of the spatial normalized ACF of intensity fluctuations [2]. Its estimate, $C (m_{x} Δ, m_{y} Δ), m_{x} = 0.. M_{x}, m_{y} = 0.. M_{y}$ , is given by:

C (m_{x}, m_{y}) \equiv C (m_{x} Δ, m_{y} Δ) = \frac{\sum_{k = 1}^{N_{x} - m_{x}} \sum_{l = 1}^{N_{y} - m_{y}} (I_{k l} - \hat{I}) (I_{k + m_{x}, l + m_{y}} - \hat{I})}{\sum_{k = 1}^{N_{x}} \sum_{l = 1}^{N_{y}} {(I_{k l} - \hat{I})}^{2}}, \hat{I} = \frac{1}{N_{x} N_{y}} \sum_{k = 1}^{N_{x}} \sum_{l = 1}^{N_{y}} I_{k l} .

This ACF was calculated only at positive spatial lags from a single speckle pattern; that’s why we omitted the index “i” for time from

I_{k l, i}

in the formulas given above. The exemplary speckle patterns, the 2D distributions of

C (m_{x} Δ, m_{y} Δ), m_{x} = 0, 1...4, m_{y} = 0, 1...4

, and the histograms of intensity for the three speckle sizes are shown in Fig. 2.

Fig. 2 From left to right: exemplary patterns for speckle types A, B, and C and 2D distributions of the spatial normalized ACF $C (m_{x} Δ, m_{y} Δ)$ of intensity fluctuations at positive spatial lags (a). Histograms of intensity for speckle types A, B and C (b).

Download Full Size | PDF

The time sequences formed at the points in Z₁ contained a lot of non-correlated entries due to $T / τ_{c}$ >> 1 whereas the sequences in Z₄ consist of practically correlated intensity values, especially for the case of $T_{1}$ . Different information capacity of the acquired statistical data in the four regions of the synthetic object affects statistical properties of the pointwise estimates given by Eqs. (1)-(4). This means that accuracy of an estimate is not the same across the activity map. It is higher in the regions of higher activity. The increase in the acquisition time improves the contrast of the built activity map. As it should be expected, the contrast of the activity map for the GD estimate is much worse than for the other three estimates. The regions Z₁, Z₂, Z₃ and Z₄ are barely recognizable for this estimate especially at $T_{1} = 64 Δ t$ and $T_{2} = 128 Δ t$ while for the other three algorithms they are clearly discernible. This is observed for the three types of speckle. For illustration, Fig. 3 depicts the activity maps obtained for the speckle “C” at $T_{3} = 256 Δ t$ . The length, $t_{w}$ , of the window used for the MGD map and the time lag $τ$ for the SF and MSF maps are equal to 4 samples. The maps were normalized to the maximum observable value in each case.

Fig. 3 Normalized activity maps obtained for different estimates for speckle C at $T_{3} = 256 Δ t$ , $t_{w} = 4 Δ t$ and $τ = 4 Δ t$ .

Download Full Size | PDF

The constructed synthetic object is very useful for comprehending the issue with the strong fluctuations of the intensity-based pointwise estimates across the activity map. We give in Fig. 4 (a) a 3D presentation of the MSF estimate across this object at $τ = 4 Δ t$ and $T_{3} = 256 Δ t$ . Large fluctuations of the estimate values due to the finite length of the temporal averaging interval T and their low spatial correlation because of processing speckle patterns are clearly seen. In view of the constant activity within the spatial regions Z₁, Z₂, Z₃ and Z₄, let us replace the MSF estimate at each point with its average value within the corresponding region. The result is depicted in Fig. 4 (b) and shows how an ideal activity map must look like, e.g. for an infinite temporal interval T. Therefore, the fluctuations of the estimate from point to point may be considered as noise.

Fig. 4 Distribution of the MSF estimate (a) and its average value (b) across the activity map of the synthetic object for the speckle type C at $T_{3} = 256 Δ t$ and $τ = 4 Δ t$ .

Download Full Size | PDF

Four histograms can be built from the entries in the four activity regions of the synthetic object. To make possible comparison of histograms for different algorithms, we transformed each 2D contour activity map to an 8-bit encoded image. Figure 5 depicts the histograms obtained when applying the four algorithms to the points in the highest activity region Z₁ for the speckle types A, B and C. The speckle size affects most strongly the GD and SF algorithms. The higher values of the estimates are observed at smaller speckle size. The histograms obtained for the regions Z₁, Z₂, Z₃ and Z₄ with the four algorithms for the speckle type C at $T_{3} = 256 Δ t$ are given in Fig. 6. Actually, these histograms correspond to the activity maps in Fig. 3. The histograms characterize behavior of the non-normalized estimates. One should bear in mind that for the normalized estimates they could be quite different. The difference arises from the pointwise normalization which uses the data from the same time sequence, $I_{k l, n} \equiv I (k Δ, l Δ, n Δ t), n = 1.. N$ , to evaluate the normalizing factor. Statistical properties of the estimates with normalization are beyond the scope of this paper, but we plan to study this issue as a future task. As it should be expected, the modes of the PDFs for the estimates take different values from 0 to 255 in the regions of different activity. The spread of the histograms for the GD estimate in Fig. 6 is very large, and the histograms in the neighbouring regions practically overlap. An overlap is observed even for distributions corresponding to the highest and the lowest activity. The histogram in the region of the lowest activity is spread almost from 0 to 255. This explains the issue with very low sensitivity of this algorithm.

Fig. 5 Histograms of estimates in the spatial region Z₁ ( $τ_{c} = 10 Δ t$ ) of the synthetic object for speckle types A (blue), B (red) and C (green) at $T_{3} = 256 Δ t$ , $t_{w} = 4 Δ t$ and $τ = 4 Δ t$ .

Download Full Size | PDF

Fig. 6 Histograms of estimates produced by GD, MGD, SF and MSF algorithms for the speckle type C in the four spatial activity regions of the synthetic object at $T_{3} = 256 Δ t$ , $t_{w} = 4 Δ t$ , $τ = 4 Δ t$ ; the numbers next to the color symbols give the correlation radius $τ_{c}$ .

Download Full Size | PDF

The results are much better for the MGD, SF and MSF algorithms. The histograms of their estimates are much narrower. It is intersting to mention that contrary to the GD estimate the spread of fluctuations of the MGD, SF and MSF estimates decreases in the regions of lower activity. This effect is most pronounced for the SF estimate. The degree of histogram narrowing for these algorithms at increase in the acquisition time is shown in Fig. 7 which presents the histograms for the lowest and highest activity regions at $T_{1} = 64 Δ t$ and $T_{4} = 512 Δ t$ . Another interesting result is that the histograms of the MGD and MSF algorithms are very close. They practically coincide for the given example. For the case of the GD estimate, the large degree of overlapping combined with a much higher computational effort excludes this algorithm from competition with the other three algorithms. That’s why we compare below only the MGD, SF and MSF algorithms.

Fig. 7 Histograms of estimates produced by GD, MGD, SF and MSF algorithms in the spatial regions with 8 times difference in the temporal correlation radius of activity at $T_{1} = 64 Δ t$ and $T_{4} = 512 Δ t$ , $t_{w} = 4 Δ t$ and $τ = 4 Δ t$ ; full symbols - $τ_{c} = 10 Δ t$ , - hollow symbols $τ_{c} = 80 Δ t$ .

Download Full Size | PDF

4. Performance assessment of activity maps from synthetic and real data

It is important to develop an approach for quantitative characterization of sensitivity of algorithms in dynamic speckle analysis. Location of the maximum in the frequency count distribution for any estimate depends on activity. Let us consider two adjacent equal size spatial regions of different activity and let us imagine a hypothetical situation when the estimate of activity exhibits the same spread of fluctuations in these regions. Such a case is depicted in Fig. 8 (a) where we have taken as an example two normal distributions to represent the PDFs of the estimate values. The areas under both curves $S_{1}$ and $S_{2}$ are equal, and the important question is how small should the area of their overlap be in order to register two differing activities. If one accepts a conventional rule that separation of the distributions requires the minimum value (equal in this case to 0.05) of one of them to coincide with the maximum value (equal to 1) of the other distribution, the overlap area $S_{12}$ should be less than 20% of $S_{1}$ .The same consideration could be applied to distributions of real estimates. For example, Fig. 8 (b) and Fig. 8 (c) show the overlap in the case of the GD and MSF estimates for the histograms obtained in the regions Z₂ and Z₃ of the synthetic object at $T_{3} = 256 Δ t$ . The almost 80% of overlap for the GD estimate makes this algorithm practically insensitive to activity change from $τ_{c} = 20 Δ t$ to $τ_{c} = 40 Δ t$ . Of course, one should keep in mind that $τ_{c}$ governs the change in activity through the exponential ACF $R (τ) = σ^{2} \exp (- τ / τ_{c})$ . The MSF performance at $τ = 4 Δ t$ is much better.

Fig. 8 Definition of sensitivity criterion for an estimate with the same normal PDF in the regions of different activity (a). Histograms of the GD estimate (b) and MSF estimate (c) in the neighboring regions Z₂ and Z₃ of the synthetic object at $T_{3} = 256 Δ t$ , $τ = 4 Δ t$ .

Download Full Size | PDF

One may argue that the requirement of 20% overlap is too strong in the case of dynamic speckle analysis, especially in view of the very low spatial correlation between the values of a given estimate within the activity map. The latter means that a relatively small spatial region in this map provides enough data to identify a given activity. The size of this region is directly related to the overlap area of the PDF curve which describes the estimate at this particular activity with the PDF curves corresponding to activity observed in the adjacent spatial regions. The smaller the overlap is, the smaller the required size is and the higher the spatial resolution of the method is as well. Therefore, the ratio

α = \frac{S_{12}}{S_{1}} \times 100 %

can be accepted as a valuable parameter for quantitative comparison of the algorithms and for evaluation of their spatial resolution and sensitivity. Note that for any two histograms built from the same number of counts at the same bin size one has

S_{1} = S_{2}

.

We calculated the parameter $α$ by numerical integration of the overlapped area for the MGD, SF and MSF estimates in the neighboring spatial regions of the synthetic object. The results are shown in Fig. 9 (a) and Fig. 9 (b) which give the decrease of $α$ with the acquisition time when the histograms in the regions $Z_{1}$ and $Z_{2}$ , $Z_{2}$ and $Z_{3}$ , $Z_{3}$ and $Z_{4}$ are compared. The presented results correspond to speckle types A and C. Simulation proves that performance of the algorithms is better for the smaller average speckle size. The SF algorithm shows the best performance in separating the regions $Z_{2}$ , $Z_{3}$ and $Z_{3}$ , $Z_{4}$ while the MGD algorithm is the best in indicating different activity in $Z_{1}$ and $Z_{2}$ . The values obtained for the parameter $α$ are so close for the three algorithms that any of them can be used for building the activity map.

Fig. 9 Parameter $α$ as a function of the number of speckle patterns for comparing the regions $Z_{1}$ and $Z_{2}$ , $Z_{2}$ and $Z_{3}$ , $Z_{3}$ and $Z_{4}$ for speckle types A (a) and C (b); full circle – MSF, hollow circle – SF, rectangle – MGD; $t_{w} = 4 Δ t$ and $τ = 4 Δ t$ .

Download Full Size | PDF

Another important issue is the optimum choice of the running window $t_{w}$ or the time lag $τ$ . In the case of the SF and MSF algorithms, the maps at larger time lags, $m Δ t$ , are obtained at lower correlation of the input data but the number of intensity values used to build the estimate also decreases as $N - m$ . In principle, one can build a set of maps of the SF estimates at increasing time lags that could provide information about the activity time scales in space. Figure 10 shows the scale of variation of the normalized SF $S = 2 {1 - \exp (- m Δ t / τ_{c})}$ in the four activity regions of the synthetic object. We see that for lags up to $20 Δ t$ , an increase in $τ$ should lead, at least theoretically, to larger separation of modes of the PDFs corresponding to the SF estimate at different activity. However, the rise in the SF estimate fluctuations with $τ$ due to lower correlation between the compared intensity values exceeds the increase in modes separation. As a result, the contrast of activity map goes down. The fluctuations are increasing up to a time lag, which guarantees lack of correlation. Above this value of $τ$ , the spread of the histograms remains more or less the same. Therefore the choice of $τ$ should be less than the minimum value of the correlation radius of intensity fluctuations. For the MGD case, the spread of fluctuations also goes up with $t_{w}$ . We studied the effect of larger values of $t_{w}$ and $τ$ by processing the speckle patterns for our synthetic object. The results are shown in Fig. 11 as plots for the parameter $α$ for comparing the histograms in the regions $Z_{1}$ and $Z_{2}$ , $Z_{2}$ and $Z_{3}$ , $Z_{3}$ and $Z_{4}$ for the MGD, SF and MSF algorithms. The MGD algorithm definitely outperforms the other two in the higher activity regions by providing the lowest $α$ value while the SF algorithm is the best for indication of lower activity regions. As a whole, the three algorithms show more or less the same performance even at the larger $t_{w}$ and $τ$ .

Fig. 10 Change of the normalized structure function with the time lag $τ = m Δ t$ at different temporal correlation radii.

Download Full Size | PDF

Fig. 11 Parameter $α$ as a function of the number of samples in the running window, $t_{w}$ , and the time lag, $τ$ , for comparing the regions $Z_{1}$ and $Z_{2}$ , $Z_{2}$ and $Z_{3}$ , $Z_{3}$ and $Z_{4}$ ; full circle – MSF, hollow circle – SF, rectangle – MGD, $T_{3} = 256 Δ t$ .

Download Full Size | PDF

The gain in using the SF and MSF algorithms is in computation time as it can be seen from Table 1 which shows the times needed to compute SF and MSF estimates, $t_{S F}$ and $t_{M S F}$ , for speckle patterns of size 256 × 256 pixels as a fraction of the time, $t_{M G D}$ , for calculation of the MGD estimate at $τ = 4 Δ t$ and $t_{w} = 4 Δ t$ . We see that using the SF and MSF algorithms ensures two times faster processing.

Table 1. SF and MSF algorithms computation time compared to the MGD computation time

View Table

We compared the quality of activity maps of the MGD, SF and MSF algorithms by applying them to sequences of speckle patterns acquired for real objects. The first of them was a specially designed circular metal object, which contained two hollow regions of the same depth – a central circular section and an annular region – and two flat annular regions. To create different activity on its surface, the test object was covered with a transparent polyester paint to form a flat layer. In this way, the circular and annular hollow regions contained larger quantity of paint than was deposited on the flat surface of the other two annular regions. The drying of the paint produced a dynamic speckle after illumination of the object with laser light. The second object was a pea leaf. We used the experimental set-up shown in Fig. 1 with a He-Ne laser. The acquisition of the speckle patterns was made at a rate of $Δ t$ = 250 ms between the frames. The size of the processed images was 500 × 780 pixels. The estimates were calculated for a series of N = 170 images.

The activity maps obtained by applying MGD, SF and MSF algorithms at $t_{w} = 4 Δ t$ or $τ = 4 Δ t$ are shown in Fig. 12. All maps were normalized to their maximum values. For the SF estimate, the level “1” in the normalized activity map corresponds to a value below the maximum value of the estimate. This mapping does not change visualization of varying activity across the object but enhances visibility of the low activity regions at the periphery of the pea leaf. The maps show a good contrast. There is practically no difference between the results for the MGD and MSF algorithms. The areas corresponding to the hollow regions of the metal object are sharply outlined and show the same activity within them with exception of some points of higher or zero activity. The higher activity values at these points are due to increased reflection. The latter causes at some points saturation of the optical sensor when all data in the time sequence are equal to 255, and all estimates at such points become equal to zero. Activity in the hollow regions was less than on the flat surface. The leaf was put on a sheet of a thick paper, which acted as non-changing background; the values of the estimates in the background region were equal to zero.

Fig. 12 2D normalized contour maps of the MGD, SF and MSF estimates that indicate activity in a metal object covered with a polyester paint and a pea leaf.

Download Full Size | PDF

To prove the efficiency of the applied algorithms to differentiate regions with different activity, we compared a randomly chosen speckle pattern captured for the metal object to the contour map of the MSF estimate. To ensure uniform illumination and reflectivity, we took only the region of 580 × 120 pixels in the upper part of the object. This part, marked with a dash line on one of the maps in Fig. 12, contains two regions of activity corresponding to a flat region and a hollow region. The speckle pattern and the MSF map for this part of the object are shown in Fig. 13 (a) and Fig. 14 (a). As it can be seen, the hollow region is practically merged with the flat one on the speckle pattern. This means that the spatial fluctuations of intensity in both regions are the same and so the regions are non-distinguishable without statistical processing. This is confirmed by the histograms of intensity values that are built for two rectangular patches with 5412 pixels each in the flat (F) and the hollow (H) regions. On the contrary, the hollow region is clearly distinguishable from the flat one on the MSF map due to the different speed of processes developing in them. Again we built the histograms of the MSF estimate for the identical patches F and H consisting of 5412 pixels and obtained two different distributions in this case.

Fig. 13 Part of a bmp speckle pattern acquired for the metal object (a), histograms of intensity fluctuations in the spatial patches F and H.

Download Full Size | PDF

Fig. 14 2D contour map of the MSF estimate for the part of the metal object (a); histograms of the MSF estimate in the spatial patches F and H.

Download Full Size | PDF

5. Conclusions

We studied performance of intensity-based pointwise algorithms for 2D evaluation of activity in optical metrology with dynamic speckle analysis. The algorithms are applied to a time sequence of correlated speckle images of a changing object illuminated with laser light. They yield a map of activity by averaging over time sequences formed from intensity values in the acquired images at the same spatial points. We proposed performance assessment of the algorithms by building histograms of their estimates. Such an approach can be applied to any pointwise algorithm, which operates with intensities or with their spectra and enables evaluation of the algorithm sensitivity to variable activity across the object. In this paper, however, we evaluated performance only of the algorithms without pointwise normalization. As a rule, they provide high contrast maps, but require uniform illumination and the same reflectivity all over the object surface. We compared characterization of activity by the GD of intensities with and without a running window to computation of a temporal SF. We showed that the window in the GD technique or the time lag for calculation of the SF narrow the spread of spatial fluctuations of the estimates within the activity map. Both approaches produce activity maps of comparable quality. Their quantitative comparison was done by analyzing the histograms of the estimates, which were built from the speckle patterns acquired for a synthetic object with four regions of different activity. As a new statistical measure, we proposed computation of a modified SF which provided the same contrast of the activity map as the other algorithms but at less computational effort. The compared algorithms were applied also to processing of speckle patterns from real objects. Processing of speckle patterns results in strong fluctuations of any estimate across the activity map and in contrast worsening due to overlap of histograms of the estimate values for the regions of different activity. Thus improvement of the contrast and sensitivity of any intensity-based algorithm can be achieved by filtering the activity map that will lead to narrowing of the spread of fluctuations. Obviously, analysis of fluctuations of the estimate used to build the activity map will facilitate the choice of a proper filter. We plan to design suitable filters as our future task.

Acknowledgment

The work is supported by the Academy of Finland, grant No. 137012: High-Resolution Digital Holography: A Modern Signal Processing Approach.

References and links

1. H. J. Rabal and R. A. Braga, Jr., eds., Dynamic Laser Speckle and Applications (CRC Press, 2009).

2. J. Goodman, Speckle Phenomena in Optics: Theory and Applications (Roberts and Company Publishers, 2007).

3. B. Zheng, C. M. Pleass, and C. S. Ih, “Feature information extraction from dynamic biospeckle,” Appl. Opt. 33(2), 231–237 (1994). [CrossRef] [PubMed]

4. R. Arizaga, N. Cap, H. Rabal, and M. Trivi, “Display of the local activity using dynamical speckle patterns,” Opt. Eng. 41(2), 287–294 (2002). [CrossRef]

5. V. Rajan, B. Varghese, T. G. van Leeuwen, and W. Steenbergen, “Speckles in laser Doppler perfusion imaging,” Opt. Lett. 31(4), 468–470 (2006). [CrossRef] [PubMed]

6. A. Serov and T. Lasser, “High-speed laser Doppler perfusion imaging using an integrating CMOS image sensor,” Opt. Express 13(17), 6416–6428 (2005). [CrossRef] [PubMed]

7. S. E. Murialdo, G. H. Sendra, L. I. Passoni, R. Arizaga, J. F. Gonzalez, H. Rabal, and M. Trivi, “Analysis of bacterial chemotactic response using dynamic laser speckle,” J. Biomed. Opt. 14(6), 064015 (2009). [CrossRef] [PubMed]

8. Z. Xu, C. Joenathan, and B. Khorana, “Temporal and spatial properties of the timevarying speckles of botanical specimens,” Opt. Eng. 34(5), 1487–1502 (1995). [CrossRef]

9. R. A. Braga, L. Dupuy, M. Pasqual, and R. R. Cardoso, “Live biospeckle laser imaging of root tissues,” Eur. Biophys. J. 38(5), 679–686 (2009). [CrossRef] [PubMed]

10. R. A. Braga, I. M. Dal Fabbro, F. M. Borem, G. Rabelo, R. Arizaga, H. J. Rabal, and M. Trivi, “Assessment of seed viability by laser speckle technique,” Biosystems Eng. 86(3), 287–294 (2003). [CrossRef]

11. R. Braga Jr, G. Rabelo, L. Granato, E. Santos, J. Machado, R. Arizaga, H. Rabal, and M. Trivi, “Detection of fungi in beans by the laser biospeckle technique,” Biosystems Eng. 91(4), 465–469 (2005). [CrossRef]

12. R. S. Macedo, J. B. Barreto Filho, R. A. Braga Jr, and G. F. Rabelo, “Sperm motility decreasing and semen fertility in the bull evaluated by biospeckle,” Reprod. Fertil. Dev. 22(1), 170–171 (2009). [CrossRef]

13. M. F. Limia, A. M. Núñez, H. Rabal, and M. Trivi, “Wavelet transform analysis of dynamic speckle patterns texture,” Appl. Opt. 41(32), 6745–6750 (2002). [CrossRef] [PubMed]

14. M. Z. Ansari, P. D. Minz, and A. K. Nirala, “Fruit quality evaluation using biospeckle techniques,” in Proceedings of IEEE Conference on Recent Advances in Information Technology (IEEE, 2012), pp. 873–876.

15. C. Mulone, N. Budini, F. M. Vincitorio, C. Freyre, A. J. López Díaz, and A. R. Rego, “Analysis of strawberry ripening by dynamic speckle measurements,” Proc. SPIE 8785, 87851X (2013). [CrossRef]

16. T. Lyubenova, E. Stoykova, E. Nacheva, B. Ivanov, I. Panchev, and V. Sainov, “Monitoring of bread cooling by statistical analysis of laser speckle patterns,” Proc. SPIE 8770, 87700S (2013). [CrossRef]

17. R. Braga, B. Silva, G. Rabelo, R. Costa, A. Enes, N. Cap, H. Rabal, R. Arizaga, M. Trivi, and G. Horgan, “Reliability of biospeckle image analysis,” Opt. Lasers Eng. 45(3), 390–395 (2007). [CrossRef]

18. J. Moreira, R. R. Cardoso, and R. A. Braga, “Quality test protocol to dynamic laser speckle analysis,” Opt. Lasers Eng. 61(17), 8–13 (2014). [CrossRef]

19. K. M. Ribeiro, R. A. B. Júnior, T. Sáfadi, and G. Horgan, “Comparison between Fourier and Wavelets transforms in biospeckle signals,” Appl. Math. 4(11), 11–22 (2013). [CrossRef]

20. A. V. Saúde, F. S. de Menezes, P. L. Freitas, G. F. Rabelo, and R. A. Braga Jr., “Alternative measures for biospeckle image analysis,” J. Opt. Soc. Am. A 29(8), 1648–1658 (2012). [CrossRef] [PubMed]

21. E. Stoykova, B. Ivanov, and T. Nikova, “Correlation-based pointwise processing of dynamic speckle patterns,” Opt. Lett. 39(1), 115–118 (2014). [CrossRef] [PubMed]

22. A. V. Saúde, F. S. de Menezes, P. L. S. Freitas, and G. F. Rabelo, “On Generalized Differences for Biospeckle Image Analysis,” in Proceedings of IEEE Conference on Graphics, Patterns and Images (IEEE, 2010), pp. 209–215.

23. H. Fujii, K. Nohira, Y. Yamamoto, H. Ikawa, and T. Ohura, “Evaluation of blood flow by laser speckle imaging sensing Part I,” Appl. Opt. 26(24), 5321–5325 (1987). [CrossRef] [PubMed]

24. A. Federico, G. H. Kaufmann, G. E. Galizzi, H. Rabal, M. Trivi, and R. Arizaga, “Simulation of dynamic speckle sequences and its application to the analysis of transient processes,” Opt. Commun. 260(2), 493–499 (2006). [CrossRef]

25. E. Equis and P. Jacquot, “Simulation of speckle complex amplitude: advocating the linear model,” Proc. SPIE 6341, 634138 (2006). [CrossRef]

N	64	128	256	512
$t_{S F} / t_{M G D}$	0.39	0.501	0.555	0.567
$t_{M S F} / t_{M G D}$	0.227	0.398	0.485	0.551

Performance of intensity-based non-normalized pointwise algorithms in dynamic speckle analysis

Abstract

1. Introduction

2. Pointwise intensity-based algorithms without normalization

3. Histogram of the estimate fluctuations

4. Performance assessment of activity maps from synthetic and real data

5. Conclusions

Acknowledgment

References and links

Cited By

Figures (14)

Tables (1)

Equations (7)

Optics Express