1.

Introduction

Laser speckle contrast imaging (LSCI) is a full-field optical technique used for imaging blood flow changes in vivo without scanning.^{1, 2, 3, 4, 5, 6, 7} The speckle pattern is generated by random interference of the scattered laser light and is modulated by motion of the scatterers. LSCI accesses the speed information of the scatterers by calculating the speckle contrast $\left(K\right)$ ,⁸ which is defined as the ratio of the standard deviation of intensity to the mean value of intensity in the speckle patterns.

As illustrated in Fig. 1 , the laser speckle contrast analysis can be performed based on spatial statistics,² temporal statistics,^{9, 10} or a combination of both.^{3, 11, 12} The laser speckle spatial contrast analysis (LSSCA) method² performs speckle contrast calculation in the spatial domain using a spatial window. LSSCA achieves high temporal resolution with the loss of spatial resolution, impeding its application on monitoring blood flow changes in small vessels. The laser speckle temporal contrast analysis (LSTCA) method,^{9, 10} which is based on temporal statistics, computes speckle contrast images using a sequence of speckle images acquired along a few time points instead of using a spatial window. LSTCA preserves the original spatial resolution by sacrificing the temporal resolution, making it inappropriate in applications where video frame rate visualization of blood flow is required. To compromise spatial and temporal resolution, several methods that attempted to combine spatial and temporal statistics have been presented. For each of these methods, a total number of $N_s\times N_s\times N$ pixels is used for calculating a $K$ value, where $N_s$ is the length of a square matrix in the spatial domain and $N$ is the number of frames. The spatial-based laser speckle contrast analysis (SLASCA) method³ is an improvement of LSSCA, in which $N$ frames of speckle contrast images computed by LSSCA with a $N_s\times N_s$ spatial window are temporally averaged. By performing a temporal average, SLASCA displays much lower image noise than LSSCA does. The temporal-based laser speckle contrast analysis (TLASCA) method¹¹ is a modification of LSTCA, in which a speckle contrast image computed by LSTCA from $N$ frames of raw speckle images is spatially averaged using a $N_s\times N_s$ square matrix (such as $3\times 3\text{pixels}$ ). By performing a spatial average, TLASCA obtains a speckle contrast image with acceptable signal-to-noise ratio (SNR) with fewer number of frames than LSTCA does. Recently, a spatiotemporal laser speckle contrast analysis (STLASCA) method¹² has been presented, in which speckle contrast is calculated directly from an $N_s\times N_s\times N$ pixel cube. Because all the $N_s\times N_s\times N$ pixel are used to compute a $K$ value, STLASCA achieves high SNR without a further averaging procedure as TLASCA or SLASCA does.

Fig. 1

Schematic diagrams of methods for laser speckle contrast analysis based on (a) spatial statistics, (b) temporal statistics, and (c) a combination of both. The $R\times C$ arrays represent the original speckle images. The elements with bright red color represent the pixels involved in the calculation of the local speckle contrast of the pixel with black color. The blue colored squares represent the derived speckle contrast images. (Color online only.)

However, the focus of previous studies has been on the processing time,^{11, 12} but the statistical accuracy of these methods, which is important in LSCI, has not yet been thoroughly compared. Since the introduction of some parallel data processing devices such as the graphics processing unit into LSCI,¹³ the computational efficiency has been significantly promoted, thus the data processing time should no longer be the primary consideration. In contrast, the close correlation of speckle contrast value with speed of motion determines that the statistical accuracy of a method for computing speckle contrast is very important. Through numerical simulation and in vivo rat cortical blood flow imaging, this study investigates the mean speckle contrast values and the relative noises of the speckle contrast images computed by SLASCA, TLASCA, and STLASCA.

2.

Methods and Materials

2.2.

Data Analysis of the Simulated Time-Integrated Dynamic Speckle Sequence

One frame of the simulated time-integrated dynamic speckle images is shown in Fig. 2 . Due to the effect of time integration, the speckle image is to some degree blurred. Figure 2 shows the intensity probability density function (PDF) of one of the simulated fully developed speckle images (green solid line), the theoretical negative exponential intensity PDF curve for a fully developed speckle image⁸ (black dashed line), and the intensity PDF of one of the simulated time-integrated dynamic speckle images (red solid line). As is shown, the intensity PDF of the simulated fully developed speckle image agrees with the theoretical negative exponential curve very well, demonstrating the validity of the simulation method. As is also shown, intensity PDF of the simulated time-integrated dynamic speckle image is no longer negative-exponential shaped due to time integration.

Fig. 2

(a) Simulated time-integrated dynamic speckle pattern. (b) Intensity probability density function (PDF) of the simulated time-integrated dynamic speckle images (red solid line), intensity PDF of the simulated fully developed speckle image (green solid line), and intensity PDF of the theoretical negative exponential intensity PDF curve for a fully developed speckle image (black dashed line). (Color online only.)

$K$ values for these simulated speckle images are then computed using SLASCA, TLASCA, and STLASCA. The $N_s$ and $N$ values are variable. Each method obtains one frame of $K$ image using a certain $N_s$ and an $N$ value. For each method, the ${\mu }_K$ and the ${\sigma }_K$ values, which are the mean value and standard deviation of $K$ , respectively, are calculated from the whole $K$ image. Subsequently, the ${\sigma }_K∕{\mu }_K$ value, which is used to quantify the noise level of a speckle contrast image, is obtained. Global speckle contrast values, which are calculated over all the pixels of the speckle image realizations, are obtained as references of speckle contrast values.

2.3.

In Vivo Rat Cortical Blood Flow Imaging

An in vivo cortical blood flow imaging experiment is also performed to compare these three methods. The schematic setup for the experiment is shown in Fig. 3 . An adult male Wistar rat weighing around $200\mathrm{g}$ was anesthetized and fixed in a stereotaxic instrument. An approximately $5.0\times 5.0\text{-}\mathrm{mm}$ cranial window with intact dura was formed by removing the skull overlying one side of the parietal cortex with a high speed dental drill (Fine Science Tools, USA) under constant saline cooling. A beam of He-Ne laser (Melles Griot, America; $632.8\mathrm{nm}$ and $15\mathrm{mW}$ ) was expanded and collimated to illuminate the cranial window at about $30\text{-}\mathrm{deg}$ incidence. 30 frames of statistically independent laser speckle images were acquired by a $12\text{-bit}$ charge-coupled device (CCD) camera (PixelFly QE, PCO Computer, Germany; pixel $\text{size}=6.45\times 6.45\mu \mathrm{m}$ ) attached to a microscope (Z16 APO, Leica, Germany; working distance $97\mathrm{mm}$ ) for data processing. The CCD exposure duration was $20\mathrm{ms}$ and the frame interval time is approximately $87\mathrm{ms}$ . The system magnification is adjusted to $3.15\times$ , and the aperture diaphragm is well controlled to ensure the average speckle size of the images to be approximately two pixels.¹⁶ A variable attenuator was used in the light path to ensure the light intensity within the dynamic range of the CCD camera. The whole setup was placed on a vibration-isolator table (VH3036W, Newport).

Fig. 3

Schematic setup for laser speckle contrast imaging.

2.4.

Data Analysis of Speckle Images Obtained from In Vivo Rat Cortical Blood Flow Imaging

The data were processed by programs written in MATLAB language. An image of the cortex under white light illumination is shown in Fig. 4 . A typical speckle contrast image computed by STLASCA over a $3\times 3\times 20\text{pixels}$ stack is shown in Fig. 4. The broken vessels that can be identified in Fig. 4 disappear in Fig. 4, demonstrating again that laser speckle contrast analysis obtains the speed information of the scatterers, rather than simply the structure pattern of the vessels.¹⁷ Two small rectangular ROIs, $R1$ with $157\times 151\text{pixels}$ and $R2$ with $107\times 211\text{pixels}$ of the speckle images, are selected for analysis. The magnified images of $R1$ and $R2$ are shown in Figs. 4 and 4, respectively. The red lines depict the center lines of the selected vessels. The locations of the vessel center lines are obtained by performing cubic polynomial fit of the skeleton of the speckle contrast images with removed spurs. The number of pixels from point $P1$ to $P2$ of vessel 1 and from $P3$ to $P4$ of vessel 2 along the vessel center lines are both 130, which is approximately $266\mu \mathrm{m}$ in distance. In such a short distance, the blood flow along the vessel center line can be considered as uniform, thus the speckle contrast values have no significant fluctuation. By setting $N_s$ to 3 and $N$ to values ranging from 3 to 30, each method obtains two groups of ${\mu }_K$ and ${\sigma }_K∕{\mu }_K$ values. The first group of ${\mu }_K$ and ${\sigma }_K∕{\mu }_K$ values are calculated from the 130 $K$ values along vessel center line 1, and the second group of ${\mu }_K$ and ${\sigma }_K∕{\mu }_K$ values are calculated from the 130 $K$ values along vessel center line 2.

Fig. 4

(a) White light image of the rat’s left parietal lobe with removed skull. (b) Speckle contrast image calculated by STLASCA method with a $3\times 3\times 20$ kernel. (c) and (d) are the magnified images corresponding to the two rectangular ROIs $R1$ and $R2$ in (b). The red lines depict the vessel center lines. (Color online only.)

3.

Results

The ${\mu }_K$ values as a function of number of frames $N$ are shown in Figs. 5 and 5 . The length of the square matrix $N_s$ is selected to be three pixels for Fig. 5 and nine pixels for Fig. 5. The ${\mu }_K$ values computed by STLASCA increase with either increasing $N$ or increasing $N_s$ . The ${\mu }_K$ values computed by TLASCA increase with increasing $N$ , but remain unchanged with increasing $N_s$ . The ${\mu }_K$ values computed by SLASCA increase with increasing $N_s$ , but remain unchanged with increasing $N$ . As is shown in either Fig. 5 or Fig. 5, the ${\mu }_K$ value computed by STLASCA is closest to the global contrast value for a given $N$ value. The global speckle contrast curve keeps almost a constant of 0.33. Given certain $N_s$ and $N$ values, the relative error between the mean speckle contrast value calculated by SLASCA, TLASCA, or STLASCA and the global speckle contrast value can be obtained. In Fig. 5, it is shown that the maximum relative error of ${\mu }_K$ computed by STLASCA is approximately 5% when $N$ is 3, while it is higher than 13% when using the other methods. In Fig. 5, the ${\mu }_K$ values computed by STLASCA almost coincide with global speckle contrast values, demonstrating that STLASCA utilizes thoroughly the number of pixels that are used for calculating one $K$ value.

Fig. 5

(a) and (b) The ${\mu }_K$ values and (c) and (d) the ${\sigma }_K∕{\mu }_K$ values for the simulated time-integrated dynamic speckle sequence as a function of $N$ . An $N_s$ value of 3 is used for (a) and (c), and an $N_s$ value of 9 is used for (b) and (d).

The relative noise of the $K$ images computed by the previously mentioned three methods can be quantified by ${\sigma }_K∕{\mu }_K$ . The ${\sigma }_K∕{\mu }_K$ values as a function of number of frames $N$ are shown in Figs. 5 and 5. $N_s$ values of three and nine pixels are selected for Figs. 5 and 5, respectively. As is shown, ${\sigma }_K∕{\mu }_K$ correlates with both the size of $N_s$ and $N$ . For a given $N_s$ value, ${\sigma }_K∕{\mu }_K$ values computed by all of the three methods decrease with increasing $N$ . On the other hand, for a given $N$ value, the ${\sigma }_K∕{\mu }_K$ values decrease with increasing $N_s$ . It is shown that the ${\sigma }_K∕{\mu }_K$ value computed by TLASCA is slightly higher than the one computed by the others when $N$ is less than 5. This is due to the fact that TLASCA is mainly based on temporal statistics, thus an insufficient number of frames will result in a high ${\sigma }_K$ value. With increasing $N$ , the difference in the ${\sigma }_K∕{\mu }_K$ values computed by the three methods becomes insignificant.

The ${\mu }_K$ values calculated from each group of 130 $K$ values along vessel center line 1 and vessel center line 2 as a function of number of frames $N$ are shown in Figs. 6 and 6 , respectively. $N_s$ value of 3 is used for computing $K$ values. As is shown, the ${\mu }_K$ values obtained from vessel 1 are lower than those obtained from vessel 2, which suggest that the flow speed of vessel 1 is faster than vessel 2. It is shown in Figs. 6 and 6 that the changes in the ${\mu }_K$ values as a function of $N$ computed by these three methods is consistent with the simulation results in Fig. 5. Likewise, the changes in the ${\sigma }_K∕{\mu }_K$ values as a function of $N$ in Figs. 6 and 6 are also consistent with the simulation results in Fig. 5. To save paper, results obtained with other $N_s$ values are not shown, in which changes in the ${\mu }_K$ and the ${\sigma }_K∕{\mu }_K$ values as functions of $N$ are similar to the simulation results.

Fig. 6

(a) and (b) The ${\mu }_K$ and (c) and (d) the ${\sigma }_K∕{\mu }_K$ values as a function of $N$ for the animal experiment. The ${\mu }_K$ and ${\sigma }_K$ values are computed from 130 $K$ values along the centerline of [(a) and (c)] vessel 1 and the center line of [(b) and (d)] vessel 2.

The ${\sigma }_K∕{\mu }_K$ values computed by STLASCA as a function of number of pixels $N_p$ used for computing one speckle contrast value are then investigated, where $N_p=N_s\times N_s\times N$ . The result is shown in Fig. 7 . The markers represent the ${\sigma }_K∕{\mu }_K$ values obtained by simulation and from vessels 1 and 2. The black solid line draws the expected tendency of ${\sigma }_K∕{\mu }_K$ values changing with $N_p$ . As is shown, the ${\sigma }_K∕{\mu }_K$ values are obtained from either simulation or experiment, and they scale with the number of pixels as $N_p^{-0.5}$ , which is consistent with previous reports where the ${\sigma }_K∕{\mu }_K$ value computed by LASCA or SLASCA as a function of $N_p$ is investigated.^{18, 19}

Fig. 7

The ${\sigma }_K∕{\mu }_K$ values as a function of number of pixels.

4.

Discussion

As shown in both Figs. 5 and 6, for the three methods SLASCA, TLASCA, and STLASCA, the ${\sigma }_K$ values correlate with both the $N_s$ and $N$ values, while the ${\mu }_K$ values show different changes with $N_s$ and $N$ . For SLASCA, the ${\mu }_K$ value is dependent on the $N_s$ value, but is independent of the $N$ value. Therefore, SLASCA results in the same ${\mu }_K$ value as the one obtained by the basic LSSCA method. For TLASCA, the ${\mu }_K$ value is dependent on the $N$ value, but is independent to the $N_s$ value. Therefore, TLASCA results in the same ${\mu }_K$ value as the one obtained by the basic LSTCA method. For STLASCA, on the contrary, the ${\mu }_K$ value correlates with both the $N_s$ and $N$ values used, so increase in either the $N_s$ or $N$ value results in higher ${\mu }_K$ value. The results shown in Figs. 5 and 6 demonstrate that STLASCA utilizes thoroughly the number of pixels $N_p$ , and thereby provide higher statistical accuracy in the laser speckle contrast value than SLASCA and TLASCA. For comparison in detail, suppose an $N_s$ value of 7 is used. According to the correlation between the ${\sigma }_K∕{\mu }_K$ value and $N_p$ , the relative noise is approximately 14% for a single speckle contrast image computed by LSSCA. If five frames of the computed speckle contrast images are averaged, then the relative noise is reduced to approximately 5%, which is more than 2.5 times lower than the relative noise of a single speckle contrast image, demonstrating the validation of SLASCA in reducing image noise. However, as is shown in our numerical simulation result, there remains approximately 5% relative error in ${\mu }_K$ for SLASCA. If TLASCA is used, with $N_s$ and $N$ values the same as those used for SLASCA, then the relative error in ${\mu }_K$ for TLASCA is approximately 7%, which is a little higher than SLASCA. If STLASCA is used, the relative error in ${\mu }_K$ can be reduced to be within 0.5%, which is much smaller than those obtained by SLASCA and TLASCA. As has been pointed out, the statistical accuracy of a speckle contrast analysis method is important in LSCI.¹⁶ The method that achieves maximized ${\mu }_K$ means that it maximizes the variation in $K$ , and thereby maximizes the variation in laser speckle contrast image. Our results demonstrate that STLASCA achieves more statistical accuracy than SLASCA and TLASCA do, given the same $N_s$ and $N$ values for these methods. In the practical applications of LSCI, statistical accuracy of speckle contrast, spatial resolution, and temporal resolution are three important measures that should be taken into consideration comprehensively. For practical application of STLASCA, a $3\times 3\times 15$ or $5\times 5\times 5\text{pixel}$ stack is recommended, resulting in approximately 0.8% relative error in ${\mu }_K$ and approximately 8% relative noise in the speckle contrast image.

5.

Conclusion

The statistical accuracy of a laser speckle contrast analysis method for blood flow imaging is important. This study quantitatively compares the statistical accuracy of the currently presented spatiotemporal laser speckle contrast analysis methods through both numerical simulation and in vivo laser speckle imaging of rat cortical blood flow. We demonstrate that STLASCA most effectively utilizes the number of pixels used for calculating one $K$ value, and thereby achieves higher statistical accuracy in the mean speckle contrast than TLASCA and SLASCA do, which is especially meaningful for LSCI being a potentially quantitative tool of estimating blood flow.