In 1981, Fercher and Briers proposed and demonstrated that single exposure laser speckle images could be used to indirectly assess the microcirculatory blood flow by analyzing the local speckle contrast.1 This technique, commonly referred to as laser speckle contrast analysis (LASCA) or laser speckle contrast imaging (LSCI), quantifies the magnitude of speckle movement by analyzing the spatial blurring that occurs as an effect of the Doppler shifts that are introduced when moving red blood cells scatter light. Since then the simplicity, usefulness, and applicability of the technique have been demonstrated in numerous medical studies.23.4.–5
An alternative way of quantifying the Doppler shifts is to study the temporal speckle fluctuations. In laser Doppler flowmetry (LDF), this is done by calculating the first moment of the Doppler power spectrum. It has been shown that this measure increases with both the speed and the amount of red blood cells.6,7 Similar studies have experimentally shown that the LASCA technique displays comparable trends.89.10.–11 However, for the LASCA technique, an additional dependency on the image exposure time is introduced.
The similarity between LASCA and LDF is further strengthened by the theory describing how the speckle contrast depends directly on the Doppler power spectrum12 or the field correlation function.13 According to LDF theory, there is also a direct dependency among the field correlation function, the intensity correlation function, the optical Doppler spectrum (i.e., the histogram of Doppler frequencies), and the Doppler power spectrum.14,15 Based on these relationships, Thompson and Andrews have demonstrated that the Doppler power spectrum can in fact be deduced from the speckle contrast if the contrast is resolved as a function of exposure time.16
In most LASCA studies, however, only a single exposure time has been used. As a single exposure LASCA image contains information from all Doppler frequencies, it will produce perfusion estimates that to some extent are sensitive to both the speed and the amount of moving red blood cells. Unlike the LDF approach, single exposure LASCA can only be used for studying the combined effect of speed and blood amount. How these two properties affect the speckle contrast is nontrivial and not yet fully explored, although some recent studies have focused on that.17,18 Different ways of interpreting how the contrast is linked to the true tissue perfusion have been presented. Still the question of how well these proposed perfusion measures reflect the actual tissue perfusion remains unanswered. Furthermore, it is known that the optical and geometrical properties of the tissue under study affect the estimated perfusion value for both LDF and LASCA. How these properties affect the LDF-perfusion is known to some extent,19,20 whereas this is just beginning to be explored for LASCA.21
The aim of this article was to show how perfusion estimates calculated using LASCA depend on changes in blood flow speed distribution, blood tissue fraction, and optical properties and how they relate to the LDF-perfusion estimate. In order to do that, a comprehensive theoretical framework showing the relationship between LASCA and LDF is presented. This includes direct analytical forms to calculate the contrast as a function of exposure time from the Doppler power spectrum and vice versa and numerical methods for simulating both photon transport in tissue and the formation of dynamic speckle patterns on a detector.
The theory section is outlined in Fig. 1. It is described how the optical Doppler spectrum, i.e., the distribution of Doppler shifts, can be calculated from a tissue model (arrow 1). This is done either by assuming only single shifts with a certain speed of the red blood cells (RBC) (Sec. 2.1) or from a multilayered skin tissue model with certain geometrical and optical properties, blood concentration, and speed distribution (Sec. 2.2). The Doppler power spectrum is calculated as the autocorrelation of the optical Doppler spectrum (arrow 2).7 It is then shown how the contrast as a function of exposure time can be analytically calculated from the Doppler power spectrum via the intensity and field correlation functions (arrows 3, Sec. 2.3). For comparison, speckle simulations are run based on the optical Doppler spectrum, and the contrast is then calculated from the speckle simulations (arrows 4, Sec. 2.4). To close the loop, it is also described how the Doppler power spectrum can be analytically calculated from the contrast (arrow 5, Sec. 2.5).
Analytic Calculation of Single Shifted Optical Doppler Spectrum
When the scattering phase function of blood is represented by a Gegenbauer kernel phase function22 with the second parameter , the histogram over single shifted Doppler frequencies RBC:s moving with speed can be expressed analytically. We have described this in an appendix to a previous paper,23 but it is given again here for the consistency of this paper.
The Doppler frequency shift that results when light is scattered an angle by an RBC moving with velocity can be expressed by24) and the probability of can be assumed to be rectangular distributed between 0 and 1 for positive frequency shifts. The probability density function for is given by
The probability density function for is given by the scattering phase function (Gegenbauer kernel) which with is
The probability density function, for Doppler shift with frequency is thus calculated by
In practice, the probability for Doppler frequencies should rather be divided into a histogram, , with frequency bins to
The single shifted optical Doppler spectrum was also simulated by calculating Eq. (1) for a large number of random values for and . The values of were chosen from a rectangular distribution, whereas the values of (for ) were chosen from Ref. 25
From Tissue Model to Optical Doppler Spectrum
We have previously described the details of how an optical Doppler spectrum can be calculated for a given speed distribution from a multilayered biooptical model.23,26 The method is based on path length distributions from Monte Carlo simulations from multilayered models with given geometries and scattering properties, where the effect of absorption is added using Beer–Lambert’s law. The optical Doppler spectrum is then calculated based on single shifted optical Doppler spectra [Eq. (6)] for each speed in the given speed distribution, which are multiple shifted with the number of shifts given by Poisson distributions that are based on the path length, the fraction of blood, and the scattering coefficient of blood.
In this study, Doppler power spectra were calculated in three layer models, with one bloodless epidermal layer containing melanin, and two dermis layers containing homogeneously distributed blood. The upper dermis layer had a thickness of 0.5 mm and the lower had a semi-infinite thickness, and the layers could contain different amounts of blood. For all models in this study, the refractive index in all layers was 1.4. The blood was modeled with an absorption coefficient of 0.58 and for deoxygenized and oxygenized blood, respectively, a scattering coefficient of , and a Gegenbauer kernel scattering phase function with and , resulting in an anisotropy factor of 0.991 and a reduced scattering coefficient of . These values are valid for blood with 43% hematocrit, a hemoglobin value of blood, and mean cell hemoglobin concentration of RBC.26 The optical properties are valid for 785 nm.
The base model used in Secs. 3.4 and 3.5 had an epidermal thickness of , a reduced scattering coefficient of all layers of , and a melanin concentration of 1% (in the epidermis layer) with absorption coefficient of . The blood tissue fraction was 1% in both dermis layers with an oxygen saturation of 80% and an average speed of . In the random models in Sec. 3.6, all those parameters were randomly chosen.
The speed distribution in the model consisted of 10 rectangular distributions, each distributed between 0 and twice the average speed of that distribution. The average speeds were linearly distributed in the logarithmic plane from a minimum average speed, i.e., , where and , resulting in a maximum average speed of , where is a scaling factor. See Fig. 2 for an example of a speed distribution with .
From Doppler Power Spectrum to Contrast
The Doppler power spectrum , i.e., the power spectral density of the temporal speckle intensity fluctuations, can be calculated for using7,14 With this cross correlation, both homodyne and heterodyne mixing of shifted and unshifted E-fields are taken into account. In Eq. (8), it is assumed that 15,27 This extension is given by
The intensity registered by a square-law imaging detector such as a CCD or CMOS sensor is described by
The simulated speckle patterns in this study have been generated using a total of 2000 Doppler shifted EM waves, each with a Doppler frequency () that have been randomly picked from a predefined distribution given by the optical Doppler spectrum. In addition to those 2000 EM waves, 0 to 4000 non-Doppler shifted EM waves were used. Each EM wave was assigned a random initial phase (uniformly distributed at ), simulating a fully developed speckle pattern. The position of the EM sources was randomly placed at a disc with a diameter 5 mm placed 50 mm away from the detector. The detector constituted 101 times 101 pixels, with a pixel pitch of . The development of each speckle pattern was simulated for 64 ms, with a time resolution of 0.02 ms, resulting in 3201 frames of speckle snapshots. The contrast for each speckle for different camera exposure times was calculated by adding consecutive snapshot frames. For each simulated exposure time , the speckle contrast was calculated as Ref. 30
From Contrast to Doppler Power Spectrum
For fully Doppler shifted speckle patterns, Thompson and Andrews16 have shown that the Doppler power spectrum can be calculated using the second derivative of the squared contrast with respect to the exposure time . Following their work, the second derivative of a partially Doppler shifted speckle pattern can be described by31
In commercial LDF systems, the perfusion value is calculated from the first moment of the Doppler power spectrum, i.e.,
In commercial LASCA systems, the perfusion value is calculated from the contrast at one single exposure time, for example,5
Single Shifted Optical Doppler Spectrum
The analytic calculations of optical Doppler spectra, Eq. (6), were verified to simulations of one billion frequencies (see end of Sec. 2.1). The comparison, for three different speeds of the red blood cells, can be seen in Fig. 3.
Comparison Between Simulated and Calculated Contrast
In order to verify the analytic expression to calculate contrast for various exposure times, the calculations were compared to contrast values calculated from speckle simulations. The optical Doppler spectra that the calculations and simulations were based on were calculated according to Eq. (6), with . The optical Doppler spectra were calculated for frequency bins distributed between and 25 kHz. From each of the resulting optical Doppler spectra, 2000 unique frequencies were randomly chosen. In addition, 0 to 4000 zero frequencies were added, giving fractions of Doppler shifted light, , of 33% to 100%. Contrast values for various exposure times were calculated based on the chosen frequencies according to Eqs. (8)–(13). For the same chosen frequencies, speckle simulations were performed for time instances 0, 0.02, 0.04, …, 64 ms. Contrast was then calculated from those speckle simulations by averaging data from consecutive time instances over a length corresponding to the required exposure time.
Contrast to Doppler Power Spectrum
In Fig. 5, a comparison between original Doppler power spectra and spectra that are first calculated to contrast [Eqs. (8)–(13)] and then recalculated to Doppler power spectra [Eq. (21)] is presented. The original spectra were for three different speeds (1, 5, and , respectively), all with 80% Doppler shifted light. The Doppler power spectra were calculated for frequency bins between and 25 kHz, resulting in exposure times between 0 and 82 ms, with steps of 0.02 ms.
In practice, it is hardly realistic to achieve images with exposure times in steps of 0.02 ms with today’s cameras. Something that could be possible to achieve is contrast from exposure times with 1 ms steps, i.e., 1, 2, 3, … ms exposure times. In order to calculate a Doppler power spectrum with frequencies above 1 kHz from contrasts with only those exposure times, the contrasts for intermediate exposure times, including exposure times below 1 ms, have to be interpolated. In Fig. 6(a), the contrast that is interpolated from the points 0, 1, … 82 ms is compared to the original, high-resolution contrast. A logarithmic time scale is used in the plot where it is apparent that there are deviations between the interpolated and the original contrasts foremost between 0 and 1 ms. When calculating the Doppler power spectrum from those two contrast curves, it is apparent that especially frequencies above 1 kHz become highly distorted [Fig. 6(b)]. In this example, the original optical Doppler spectrum was based on single shifts with RBC speed of , 80% Doppler shifted light. A shape-preserving piecewise cubic interpolation (MATLAB) was used. The resulting perfusion estimates [Eq. (22)], which were calculated from the Doppler power spectra in Fig. 6(b), were 242 and 527, respectively.
Blood Perfusion Effect on Contrast
It was investigated how speed distribution and blood tissue fraction affect the contrast and the LASCA-perfusion estimate calculated with Eq. (24). In order to do that, optical Doppler spectra were calculated from a skin model as described in Sec. 2.2. The effect of changed blood tissue fraction is visualized in Fig. 7, whereas the effect of blood flow speed is visualized in Fig. 8. The LASCA-perfusion is normalized with the average perfusion for each exposure time in each plot.
Effect of Geometrical and Optical Properties
The size of each single Doppler shift is directly dependent on the scattering angle, and thus the scattering phase function of blood.7 It is therefore expected that the Doppler power spectrum and the contrast are affected by changing that parameter. The effect on the perfusion estimates when changing the anisotropy factor from the original value is summarized in Table 1. This was done both for a constant scattering coefficient so that the reduced scattering coefficient was affected, and by also changing so that was kept constant at .
Effect on LDF- and LASCA-perfusion estimates when changing scattering properties of blood, normalized to first line.
|Blood scattering properties||Normalized perfusion estimates [-]|
|g (-)||μs (mm−1)||μs′ (mm−1)||LDF-perf||1/K−1||1/K2−1|
When changing the reduced scattering coefficient of the static medium, the perfusion reached a maximum for values of about . The relative change was highest for the longest exposure time. The effect is summarized in Table 2.
Effect of changed scattering properties of the static medium. Perfusion estimates are normalized to the level at μs′=2.0 mm−1.
Changing the thickness of the epidermis layer has an effect on the fraction of Doppler shifted photons, which in turn affect the contrast and perfusion estimates. The effect is summarized in Table 3.
Effect of epidermis thickness on fraction of Doppler shifted light (fD) and perfusion estimates. Perfusion estimates are normalized to thickness of 0 mm.
|Thickness (mm)||fD (-)||LDF-perf||1/K−1||1/K2−1|
Finally, the melanin content in the skin also has a small effect on the perfusion estimates. The perfusion decreased with less than 15% when changing the melanin content from 0% to 40% except for 40% melanin concentration in -thick epidermis layer. Such high melanin concentration in such thick epidermis is hardly realistic. The effect was much higher for the LASCA-perfusion estimates and higher for thicker epidermis layer (higher total amount of melanin). The results are summarized in Table 4.
Effect of melanin concentration on perfusion estimates for two different epidermis thicknesses (tEpi). Perfusion estimates are normalized to the melanin concentration of 0% for each of the two thicknesses.
|tEpi (μm)||Melanin conc. (%)||LDF-perf||1/K−1||1/K2−1|
Based on the skin model described in Sec. 2.2, 2000 random models with varied blood concentration, speed distribution, optical properties of static media, and geometries were generated. The true perfusion in the models, i.e., the actual RBC tissue fraction times their average speed, was calculated resulting in the perfusion unit [% ]. The optical Doppler spectrum and Doppler power spectrum were calculated, and the LDF-perfusion was calculated from the power spectrum [Eq. (22)]. In addition, the contrast was calculated from the power spectrum and the LASCA-perfusion was calculated from the contrast, based on both Eqs. (23) and (24). In Fig. 9, the LDF- and LASCA-perfusions are plotted against the true perfusion. The LDF- and LASCA-perfusions are normalized at 0.05% . Data were sorted with respect to the true perfusion and divided into groups with 100 models in each group, from which average values and standard deviations were calculated and plotted.
It was also studied how the perfusion estimates were affected by isolated changes in average speed or RBC tissue fraction. In order to do that, the same 2000 models as above were used, and in each of the models either the average speed was increased by 10% or the RBC tissue fraction was increased by 10%. Both these changes would ideally generate a 10% increase in the resulting perfusion estimates. The actual response in the three calculated perfusion estimates is shown in Fig. 10. Data were sorted with respect to either average speed (a) or RBC tissue fraction (b) and divided into groups with 100 models in each group, from which average values and standard deviations were calculated. The average (standard deviation) estimated perfusion increase was 9.7 (0.2)%, 5.9 (0.9)%, and 8.1 (0.9)%, respectively, for LDF-perfusion, and the LASCA-perfusions calculated with [Eq. (23)] and [Eq. (24)] when average speed was increased 10%. Corresponding numbers when RBC tissue fraction was increased 10% were 6.6 (1.2)%, 7.3 (1.6)%, and 10 (1.5)%. The average standard deviation for the groups shown in Fig. 10 was 0.12%, 0.58%, and 0.85%, respectively, for the three perfusion estimates when average speed was increased 10%. Corresponding numbers when RBC tissue fraction was increased 10% were 0.36%, 0.62%, and 0.92%.
We have described how the contrast can be calculated from the optical Doppler spectrum, via the Doppler power spectrum, and vice versa. Together with the generation of optical Doppler spectra from relevant tissue models, this can be used as a toolbox for further studies on the LASCA technique for a deeper understanding and, potentially, development of better blood perfusion measures.
It is shown in Fig. 3 that the single shifted optical Doppler spectrum can be calculated analytically for a given speed distribution and scattering phase function in an accurate and effective manner.
The results in Fig. 4 clearly show that the analytically calculated contrast values agree with the simulated contrast values. The small differences seen, with no particular trends, most likely originate from the fact that the simulations were based on a limited number of light sources randomly distributed over an area. Increasing the number of random points decreases the differences, but as these simulations are rather time consuming, they were limited to a maximum of 6000 points in this study. It should also be noted that it is important to have small time steps in the speckle simulations for correct results. For example, for correct contrast values at 0.5 ms exposure time, the time step had to be or smaller.
When calculating the Doppler power spectrum from the contrast as in Fig. 5, the first- and second-order derivatives of are numerically estimated. For short exposure times, that estimation is not exact, which explains the small difference seen in Fig. 5 for especially the highest speed. The effect is pushed toward higher frequencies when increasing the sampling rate, i.e., decreasing the time steps in the contrast curve. Although we show that it is theoretically possible to calculate the Doppler power spectrum from the contrast curve, it is probably not the best way to go in order to calculate the best possible perfusion estimate. The results in Fig. 6 show that the Doppler power spectrum becomes highly distorted when calculating it from a contrast curve containing interpolated new contrast values from the original contrast values taken at realistic multiple exposure times. Maybe the interpolation can be improved in order to keep a continuous first- and second-order derivative, thus avoiding the overtones in the Doppler power spectrum, but even then, it is a computationally demanding way to go.
Studying the results in Figs. 7 and 8 gives some fundamental insights about the effect on contrast when changing blood flow parameters in a typical case. First, those figures reveal that the contrast and thus the LASCA-perfusion is sensitive to blood concentrations and flow speeds over a wide range. Another interesting finding is that the contrast is sensitive to the shape of the speed distribution, not only to the average speed, and nonlinearities arise at different speeds depending on that shape. This is in line with the findings by Ramirez-San-Juan et al.33 Furthermore, the nonlinearities of perfusion in respect to speed depend on the blood concentration, as seen when comparing Figs. 8(c) and 8(d), where a lower blood concentration results in more severe nonlinearities in that particular example. The results are in line with a recent publication by Khaksari and Kirkpatrick,18 where the LASCA-perfusion is studied in more detail in terms of diffusive and advective flux.
We have compared two common alternative equations for calculating a perfusion estimate based on a single exposure time, Eqs. (23) and (24). In Fig. 9, it can be seen that the latter has a more linear relationship with the true perfusion, whereas the former resembles the LDF-perfusion more closely. In most of the presented results, an exposure time of 6 ms has been used. We have observed that most trends, for example, the effect of changing optical properties, are similar for other exposure times as well, but the effects are generally smaller for shorter exposure times. Furthermore, when changing optical and geometrical properties in the model, the LASCA-perfusion estimates are more affected than the LDF-perfusion estimate, see, for example, Tables 1–4.
The results in Fig. 10 are interesting and important since these techniques are often used to study relative changes in perfusion spatially and/or temporally. It can be seen that the LDF-perfusion estimate responds better to changes in average speed than either LASCA estimate. The increase in the estimated LDF-perfusion responds almost perfectly to the increase in speed in the models. This is not surprising as the LDF-perfusion, in contrast to the LASCA-perfusion, scales linearly to changes in Doppler frequencies, which directly depend on the speed of the moving scatterers. The response in the LASCA-perfusion estimates on the other hand is generally too low and is somewhat dependent on the level of the average speed (i.e., nonlinear) and is less predictable (much higher standard deviations). For changes in RBC tissue fraction, it can be seen that all perfusion estimates are nonlinear, i.e., the response is generally smaller for models with higher RBC tissue fraction. As for an increased average speed, the response for an increased RBC tissue fraction is less predictable with the LASCA-perfusion estimates. In all cases, the response for the perfusion estimate calculated with Eq. (24), i.e., , is stronger than the response calculated with Eq. (23), i.e., , which was expected.
This article has shown that there is a direct relationship between the Doppler power spectrum and the speckle contrast as a function of exposure time, utilized in LDF and LASCA, respectively. Due to practical limitations, commercially available LASCA systems use only a single exposure time when calculating the perfusion estimate. This data reduction results in a less predictable perfusion estimate than for LDF. In order to improve the perfusion estimate given by the LASCA technique at least to the level of the conventional LDF-perfusion estimate, a multiple exposure time LASCA has to be employed.34,35
This study was financially supported by the Swedish Research Council (grant no. 2014-6141) and by the CENIIT research organization within Linköping University (project id. 11.02).
Ingemar Fredriksson is an adjunct lecturer in the Department of Biomedical Engineering, Linköping University, Sweden, and an R&D optics designer at Perimed AB, Järfälla-Stockholm, Sweden. His research focuses on modeling and model-based analysis of laser speckle based and spectroscopic techniques with applications within monitoring and imaging of microcircular blood flow and metabolic processes.