1. Introduction

Laser Speckle Contrast Imaging (LSCI) is a blood flow imaging technology that is broadly applied to assess microcirculation in the brain [1,2], retina [3,4], kidney [5,6] and other organs [7]. Blood particles scatter the laser light, forming a speckle pattern recorded on the camera sensor [8–10]. When particles move, the associated speckle pattern fluctuates and, provided a sufficiently long exposure time, appears blurred when recorded by the camera. The degree of blur is quantified by calculating the image’s contrast K. Such a contrast can then be converted to the blood flow index, which characterizes particle dynamics. Contrast K is calculated as the ratio of standard deviation ($\sigma$) over mean ($\langle I \rangle$) [9,8]:

Typically the calculation is performed either over the spatial distribution of intensities in a window of 5x5 or 7x7 pixels [9] or temporal distribution of intensities over 25 consecutive frames [11], although anisotropic [12] and spatiotemporal [13] approaches were introduced as well. Both spatial and temporal contrast analyses are used to characterize the microcirculation condition. Calculating the spatial contrast provides the best temporal resolution of blood flow changes [9,14]. In contrast, temporal contrast provides the best spatial resolution to analyze structural features of vessels, for example, the diameter [15]. Some critical parameters, such as characteristics associated with the cardiac pulse wave propagation, require maximizing both resolutions and thus are unavailable or with limitations for conventional LSCI analysis. Assessing these characteristics, e.g., pulsatility and resistivity indexes and pulsatility-associated diameter changes of microvasculature, with LSCI, would allow additional insights into the downstream resistance and overall microcirculatory network and tissue conditions [3,16–19]. Furthermore it would have a significant translational value since these characteristics are often assessed clinically in larger vessels using ultrasound blood flow imaging [18].

In the present study, we propose a new method for the lossless contrast analysis of periodic oscillatory blood flow changes such as cardiac pulsatility. The method is based on calculating temporal contrast over nonconsecutive frames to preserve both maximum temporal and maximum spatial resolutions in a composed blood flow oscillation period. To validate the method and compare it to standard contrast analyses, we have used a dynamic speckle simulation and in-vivo imaging. We show that the lossless method retains maximum spatial and temporal resolutions and thus allows more accurate characterization of microcirculation features associated with cardiac pulsatility.

2. Methods

2.1 Lossless temporal contrast analysis

The proposed method is based on calculating temporal contrast over frames that belong to the same phase of the periodic signal rather than over consecutive frames. Lossless temporal contrast calculation relies on three conditions:

A priori knowledge of each frame position in the period (phase) is required to perform the lossless temporal contrast calculation. To know the position of the frame in the period and to assess the second condition, it is necessary to estimate the temporal change of the signal. Such a change can be calculated by obtaining the time course using the spatial contrast analysis averaged over the whole field of view (or a sufficiently large region of interest). Then, the time course is used to bin the raw frames from each period according to their corresponding phase, while periods with an abnormal dynamic change (e.g., artifacts due to the mouse sporadically moving) are excluded. In the present study, to automatically exclude periods affected by motion artifacts or sudden changes in the blood flow, we have introduced several rejection criteria based on the variation in the cycle features. Specifically, we have measured average intensity, average blood flow index, pulse duration, pulse magnitude, and the foot-to-foot difference in blood flow index for every individual pulse. For each of the features, we then identified the median and the standard deviation. Cycles with features deviating from the median more than two times the standard deviation were rejected. The lossless temporal contrast (ltLSCI) is then evaluated using Eq. 1) for each pixel, where $\sigma$ and $\langle I \rangle$ are calculated over the frames in the same bin (belonging to the same phase of the cycle). The analysis results in a set of contrast frames representing the dynamic change over the composed period while preserving spatial and temporal resolutions. To compare our method with the standard ones, spatial and temporal contrasts were calculated using Eq. (1). Spatial contrast (sLSCI) was calculated with a sliding window of 5 by 5 pixels. Temporal contrast (tLSCI) was calculated for each pixel in a target frame using the corresponding pixel values over 25 consecutive frames. The processing was implemented in MATLAB and is available for download [20].

The interpretation of the contrast calculated with ltLSCI will follow the same principles as conventional temporal laser speckle contrast analysis in relation to dynamics regime and static scattering effects. Since blood flow dynamics, rather than contrast, are the main interest in most studies, we converted contrast to blood flow index (BFI) using the following equation [3–5]:

2.2 Simulation

We simulated speckle dynamics from particles moving with a periodic acceleration to validate the proposed method in a noise-free (except speckle noise) environment. The model used for the simulation is an extension of the model published in [21,22]. Briefly, N = 10000 particles were evenly distributed in a volume of 10 by 10 by 0.1 millimeters (see Fig. 1(A)). Each particle is considered a source of an electric field, which superimpose on the camera sensor at a distance from the volume that generates a speckle size of around 2. The intensity of the superimposed electric field (E) in each pixel of the camera is calculated with Eq. 3):

 

Fig. 1. Simulation of the speckle dynamics (A-C) and contrast calculation (D, E). A, projection of N=10000 particles distributed in a 3D volume. The regions with a blue background represent vessels of different calibers, with particles moving linearly along the X-axis and accelerating following the pulse signal. The other regions simulate parenchyma, with particles moving according to Brownian motion and accelerating following a sinusoidal signal. B, the intensity of the simulated speckle pattern that corresponds to T=0. C, the intensity of the simulated speckle pattern integrated over T=5 ms. D, schematic representation of how spatial contrast is calculated (sLSCI) over a 5x5 pixels neighborhood. E, schematic representation of how temporal (tLSCI) and lossless temporal (ltLSCI) contrast analyses are calculated. tLSCI is calculated over 25 frames, including frames from more than one pulse. ltLSCI uses the frames that belong to the same phase over multiple pulses (the frames used to calculate the first ltLSCI frame are indicated with green lines).

Download Full Size | PDF

To simulate structural features such as blood vessels, we separate the volume into different compartments along the Y-axis (see Fig. 1(A)). The compartments alternate between two types: vessel (ordered motion, average velocity= 2.43 mm/s to simulate blood velocity in small vessels [24]) and parenchyma (unordered motion, average velocity = 0.79 mm/s, [24]). The parenchyma compartments are set to be 7 pixels wide, so the spatial contrast windows would not overlap between different compartments. The vessels’ compartments change in size from 2 to 55 pixels to simulate different magnifications. Each compartment has a number of particles proportional to the size of the compartment (the density of particles is constant), and particles cannot travel between compartments. Along the X-axis, the sensor is set to be 100 pixels. To test temporal resolution, we modulated the speed of the particles following a typical pulse waveform (see Fig. 1(E)). Roughly, such modulation was created by adding two 10 Hz sawtooth signals dephased 20 degrees and then smoothing them with a median filter. In the parenchymal compartments, the particles moved following a sine wave of 7 Hz. To compute the lossless temporal contrast over the same number of frames as conventional temporal contrast, we generated 25 periods of particles moving, as we have described above.

2.3 In vivo imaging

All experimental protocols were approved by the Danish National Animal Experiments Inspectorate and were conducted according to the guidelines from Directive 2010/63/EU of the European Parliament on the protection of animals used for scientific purposes. We used n = 2 C57Bl6 mouse (Janvier, Denmark weight, 12 weeks old). One animal was used in the acute, through-the-skull imaging experiment, where an optically transparent cranial window was installed. Before starting surgical procedures or imaging sessions, the animals were anesthetized in a chamber with 3% isoflurane mixed with oxygen at a flow of 1L/min. During the surgery and imaging sessions, isoflurane concentration was reduced to a final concentration of 1-2%. The mouse was placed on a servo-controlled heating table, maintaining their body temperature at 37$^{\circ }$C. The surgical preparation for making an optically-transparent chronic cranial window is described in detail in [25]. Briefly, an optically-transparent 4mm round glass was installed in the mouse’s skull during craniotomy surgery. The surgery was performed in the area of the barrel cortex. After surgery, the mouse was allowed to recover for 7 days before imaging. For the acute through-the-skull imaging in the second mouse, the skin was removed from the skull surface, which was then covered with 2% agarose solution and fixed with a cover glass (to reduce the amount of the reflections and keep the skull surface moisturized).

During the chronic cranial window imaging session, a holographic volume grating stabilized laser diode coupled to a polarization-maintaining fiber (785 nm, Thorlabs FPV785P) controlled with a laser driver (Thorlabs LDC210C) and temperature controller (Thorlabs TED200C) was used to deliver coherent laser light on the window. The light was shone and collected coaxially using a polarizing cube and a 5x Mitotyou objective. A field of view of 2x1mm (1024x512 pixels) was recorded with a CMOS camera (Basler aca2040-90umnir) for 10 minutes with a frame rate of 194 frames per second. The camera and its settings were the same for the through-the-skull imaging experiment, but a different lens (VZM 450i, Edmund Optics) and laser (LP785-SAV50) were used in a side-illumination setup.

3. Results

3.1 Simulation

Simulation results showing how the different contrast calculation methods affect the spatial and temporal resolutions are shown in Fig. 2). As expected, tLSCI allows better defining the edges of narrow vessels than sLSCI, which naturally blurs the boundaries of the vessel as it uses a spatial kernel for its calculation (see Fig. 2(A, B, D, E)). The sLSCI calculation captures the temporal dynamics of the signal better than the tLSCI calculation, which smooths the signal as it uses a temporal kernel (see Fig. 2(G, H)). Finally, the ltLSCI analysis preserves both the spatial and the temporal resolutions (see Fig. 2(C, F, I)). Potential improvements in temporal dynamics between ltLSCI and sLSCI for the small vessel (Fig. 2(G, I dashed line)) are associated with the fact that spatial contrast analysis mixes signals from the inside and the outside of the vessel. The 1-pixel-wide vessel is hardly visible when using sLSCI analysis but also appears blurred with tLSCI and ltLSCI analyses. This blur effect with tLSCI and ltLSCI can be explained by the simulated speckle size (which defines optical resolution) being approximately two times larger than the pixel size.

 

Fig. 2. Contrast analysis methods comparison using simulated data. A-C, contrast images obtained by the spatial (sLSCI), temporal (tLSCI), and lossless temporal (ltLSCI) contrast analysis methods, respectively. The color bar corresponds to the contrast intensity normalized for all three contrast maps. D-F, spatial profiles of images shown in A-C, which are calculated by averaging all the rows in the corresponding images. G-I, normalized ground truth (red line) and blood flow index profiles for 15 (solid blue line) and 4 (dashed blue line) pixels-wide vessels. Normalized blood flow index profiles were calculated by averaging the pixels belonging to the respective values in space and then normalizing them to the 0-1 range. Minor note: the phase shift in the temporal contrast analysis is caused by the direction of the contrast calculation window and can be compensated by shifting the signal by half of the window size (12 frames) forward.

Download Full Size | PDF

Aside from maintaining spatial and temporal resolutions, the ltLSCI method is less affected by contrast offsets due to static scattering and temporal smoothing and, therefore, should provide more accurate measurements of the relative magnitude of the flow change. sLSCI analysis is strongly affected by static scattering and speckle statistics bias [26] (typically 12 speckles over the spatial window), most commonly leading to underestimating the contrast. tLSCI is less affected by static scattering and has inherently better statistics (typically >200 decorrelation periods over the temporal window), but its ability to accurately measure and detect fast changes is limited because of the temporal smoothing effect. For example, for the large vessel shown in Fig. 2, the mean BFI measured is 22.98, 24.56, and 16.94, and the relative magnitude is 0.51, 0.078, 0.48 for ltLSCI, tLSCI and sLSCI, respectively. Compared to the ground truth signal, this produces a relative magnitude estimation error of 35%, 79%, and 27%. The error is noticeably smaller for ltLSCI and sLSCI than for tLSCI, which is explained by the temporal smoothing in the latter. The nevertheless large error in both sLSCI and ltLSCI methods can be attributed to using a simplified speckle simulation model and, therefore, insufficient speckle statistics. Insufficient speckle statistics can also explain the reason for the sLSCI BFI being noticeably different from the other two methods.

3.2 In vivo

The in-vivo experimental data analysis is consistent with the simulation results. Figure 3 shows averaged contrast images, spatial profiles, and temporal dynamics for different contrast calculation methods. The tlLSCI approach provides a similar ability to capture temporal dynamics as sLSCI analysis (see Fig. 3(E, G)) and resolves spatial features as well as tLSCI ( see Fig. 3(D, F)). The difference becomes particularly noticeable for the small ($\approx$2 pixels wide) vessel. The relative magnitude of oscillations measured with ltLSCI, tLSCI, and sLSCI are 0.246, 0.041, and 0.228 for the large vessel and 0,225, 0,037, and 0.179 for the small vessel, respectively. One can notice that the relative magnitude measured with ltLSCI and sLSCI is very similar for the large vessel (7% relative difference) but is more different for the small one (20% relative difference). Similar to the simulation example, it can be explained by the spatial contrast analysis "mixing" parenchymal and vessel signals.

 

Fig. 3. Performance comparison of contrast analysis methods in obtaining the average blood flow pulse from in-vivo data. We have selected regions V1 and V2 as examples of a narrow vessel ($\approx$2 pixels wide) and a large vessel ($\approx$10 pixels wide). A-C, contrast maps calculated with the spatial (sLSCI), temporal (tLSCI), and proposed lossless temporal (ltLSCI) contrast approaches, respectively. The color bar represents the intensity magnitude of the contrast in the images. D and F represent the vessel’s profile after averaging all the rows of the V2 regions and all the columns of the V1 regions for the tLSCI and ltLSCI methods and taking the central row in V1 or the central column in V2 for the sLSCI method. E and G represent the blood flow index considering only the pixels inside of the vessel. For the V1 region, it corresponds to 2 pixels in width times 5 pixels in length. For the V2 region, we considered 10 pixels in width and 5 pixels along the vessel. When calculating sLSCI BFI, only one line along the vessel was used for averaging as sLSCI contrast values already represent neighborhoods of 5x5 pixels. Note that the ltLSCI contrast map (C) might appear truncated due to the same color scale used for all contrast maps. This was done to highlight the difference in contrast values and does not mean that the original data is truncated.

Download Full Size | PDF

An important application for the ltLSCI approach would be measuring properties of fast periodic signals in the presence of strong static scattering, e.g., when imaging a rodent brain through the skull. Figure 4 shows an example of contrast and relative oscillation magnitude maps calculated using the three contrast analysis methods. Unlike imaging results for the optically transparent cranial window (Fig. 3), the spatial contrast appears to be noticeably more blurred and distorted due to the static scattering (see Fig. 4(A)). On the other hand, tLSCI and ltLSCI provide sharper contrast maps of similar values and quality (see Fig. 4(B, C)). Finally, neither sLSCI nor tLSCI can accurately capture relative oscillations magnitude. In the case of the conventional temporal contrast analysis, due to the temporal smoothing of the signal, the relative magnitude values decreased $\approx$10 times ($\approx$0.03 in parenchyma), and the spatial features were mostly lost. For sLSCI, spatial features of the relative magnitude appear to be more blurred, and the values ($\approx$0.25 in parenchyma) are lower than compared to ltLSCI, which can be explained by the contribution of static scattering and speckle statistics bias. The proposed ltLSCI calculation performs better than the other two, with more spatial features (the reduced relative magnitude in veins became particularly clear) and average parenchymal contrast values around 0.35.

 

Fig. 4. Through-the-skull imaging of blood flow dynamics using different contrast analysis methods. A-C, contrast maps calculated with sLSCI, tLSCI, and ltLSCI, respectively. D-F, corresponding relative magnitude maps. The relative magnitude of the contrast was calculated as the difference between the pixel’s maximum and minimum intensity values in time, normalized by its mean value.

Download Full Size | PDF

4. Conclusion

We have presented a new method that improves the analysis of laser speckle data from periodic signals. We benchmarked our method with conventional spatial and temporal contrast calculations. The proposed method better retains spatial and temporal resolutions when analyzing periodic signals and, therefore, allows more accurate analysis of topology and blood flow dynamics in small vessels. Furthermore, we have proposed that the lossless temporal contrast analysis approach can become a particularly useful tool when studying fast blood flow dynamics in media with strong static scattering (e.g., through-the-skull imaging). In such a situation, the ltLSCI method not only retains the best features of both spatial and temporal contrast analysis but also provides improved accuracy of relative magnitude measurements. Finally, we expect that the method will also benefit the study of vasoreactivity (e.g., myogenic oscillations, vasomotion, flowmotion, and the relation of the last two) and neurovascular coupling. In the latter case, we expect that the proposed method will allow a better calculation of temporal dynamics (e.g., time-to-peak) in small vessels compared to conventional analyses. While the primary purpose of the method is to calculate an average cycle, it might be modified to be used in applications where dynamic changes in the cycle features, such as heart rate variability, are of interest. Specifically, one could envision a "sliding window" version of the proposed method or a comparison of the cycles binned according to their duration rather than averaging. Without such modifications, the method should not be applied in studies where cycle-to-cycle variability is important. It is essential to note that the ltLSCI method might be more susceptible to motion artifacts or aperiodicity in the laser speckle signal. Therefore, it requires more complex data pre-processing than traditional methods and might require data rejection (as in our approach [20]) or registration, which might not be necessary for conventional approaches.