Open Access
Add to BibSonomy
Abstract
Laser speckle contrast imaging (LSCI) has gained broad appeal as a technique to monitor tissue dynamics (broadly defined to include blood flow dynamics), in part because of its remarkable simplicity. When laser light is backscattered from a tissue, it produces speckle patterns that vary in time. A measure of the speckle field decorrelation time provides information about the tissue dynamics. In conventional LSCI, this measure requires numerical fitting to a specific theoretical model for the field decorrelation. However, this model may not be known a priori, or it may vary over the image field of view. We describe a method to reconstruct the speckle field decorrelation time that is completely model free, provided that the measured speckle dynamics are ergodic. We also extend our approach to allow for the possibility of non-ergodic measurements caused by the presence of a background static speckle field. In both ergodic and non-ergodic cases, our approach accurately retrieves the correlation time without any recourse to numerical fitting and is largely independent of camera exposure time. We apply our method to tissue phantom and in-vivo mouse brain imaging. Our aim is to facilitate and add robustness to LSCI processing methods for potential clinical or pre-clinical applications.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
In the standard implementation of laser speckle contrast imaging (LSCI), an expanded laser beam is launched into a tissue and the resultant backscattered speckle is imaged by a camera [1–6]. Dynamics within the tissue cause the speckle to decorrelate over time, and a key premise of LCSI is that this decorrelation time can be inferred from measurements of the speckle contrast, either spatial [1,7] or temporal [8,9]. However, the accurate association of speckle contrast to decorrelation time is confounded by several factors. First, imaging with a camera always involves some kind of filtering, both in space (because of sensor pixelization) and in time (because of finite exposure time). Such filtering must be corrected for [10–13]. Similarly, biases arising from shot-noise and camera noise, must also be corrected for [14,15], all of which is straightforward.
Another more fundamental confounding factor comes from the limited statistics provided by any given measurement. For example, in spatial LCSI, the finite-sized windows used for spatial contrast evaluation encompass only limited numbers of speckle grains, causing an underestimate in the contrast evaluation [13,16–18]. Similarly, in temporal LSCI the finite-duration measurements may not be long enough to encompass the full ensemble statistics of the tissue. A disparity between time average and ensemble average leads to the so-called non-ergodicity problem [19,20], which also leads to an error in the temporal contrast evaluation. Corrections to the latter typically involve introducing a static component to the speckle field dynamics [6,21–23]. This static component may be real, arising for example from laser reflections from a static interface such as an intact skull [24,25] or from an external local oscillator [26]. Alternatively, it may be more phenomenological in nature and adjusted to account for the balance of the ensemble statistics not captured by the finite-duration measurements [23], in which case the strength of the static component depends on camera exposure time and must be calibrated accordingly.
Still another confounding factor that hampers the association of speckle contrast with decorrelation time is that this association is not direct, or at least it has not been in previous reports. To our knowledge, this association has always involved the numerical fitting of data to a specific theoretical model that must be decided upon in advance [21]. The most representative such model is where the temporal speckle field correlation decays as a negative exponential, such that the logarithm of the correlation scales linearly with elapsed time [1]. This model is most adapted to regimes of single-scattering unordered (SU) motion or multiple-scattering ordered (MO) motion [27]. In another model the logarithm of the correlation scales quadratically with elapsed time, which is more appropriate for single-scattering ordered (SO) motion [27–29]. In yet another model the correlation logarithm scales as the square root of elapsed time [27,30], which has been shown to be most appropriate for multiple-scattering unordered (MU) regimes where the blood vessels are densely packed and unresolvable, as in parenchymal tissue regions [31,32]. Still other models have been described in Ref. [13]. Difficulties arise when the theoretical model is not known in advance or is varying across the field of view. To address these difficulties, a dynamic light scattering Monte Carlo (DLS-MC) method [33] can be used to reveal different behaviors of the speckle field correlation decay at short and long time scales [34]. An ultra-fast camera can also be used to resolve the speckle correlation decay at a fine enough temporal resolution to perform multi-parameter fitting. This strategy called dynamic light scattering imaging was adopted in [32], making use of frame rates as high as 20kHz. To achieve such high frame rates required the use of a specialized camera and imaging could be performed over only limited fields of view, requiring careful correction for noise. Nevertheless, this technique clearly revealed the occurrence of different dynamic regimes obeying different correlation models. A more broad characterization of these correlations that can be generalized to arbitrary models thus seems desirable.
In this paper, we present an alternative method of processing temporal LSCI data that can be implemented with a conventional camera operating at much slower frames rates, enabling large field-of-view imaging. We introduce a generalized definition of the speckle correlation time that, in cases of ergodic measurements, enables a direct association of speckle contrast with speckle field dynamics that is completely model-free and does not require numerical fitting of data. Interestingly, this generalized definition is routinely used when characterizing temporal fluctuations of light [35–37] but does not seem to have been exploited yet by the LSCI community. In cases of non-ergodic measurements, that is, when background speckle dynamics are present that are so slow as to be considered static, our method is no longer completely model free. Nevertheless, it continues to enable a direct association of speckle contrast with speckle field dynamics that remains accurate and robust while requiring no numerical fitting. Our alternative LSCI data processing strategy is intended to facilitate the practical implementation of LSCI, thus reinforcing the immense power and promise of this technique.
2. Theory
2.1 Generalized correlation time for ergodic measurements
In temporal LSCI, speckle patterns emanating from a tissue sample are imaged by a camera and recorded sequentially over time. In the idealized case where the camera response time is infinitely fast, the normalized temporal intensity correlation function at any pixel is given by
where $I(t)$ is the instantaneous intensity, and $\langle..\rangle$ corresponds to an ensemble average. Unless otherwise stated, we assume that the full ensemble statistics are captured within the overall measurement time, meaning that the system may be considered ergodic and $\langle..\rangle$ may be replaced by a temporal average. In this case, $g_2(\tau )$ is related to the normalized field correlation function $g_1(\tau )$ through the well-known Siegert relation [38] Here, $\beta$ is a correction factor that accounts for spatial filtering only, due to the finite size of the camera pixels. From Eq. (2), the intensity correlation time is directly related to the field correlation time, and indeed, a measurement of the latter is one of the aims of LSCI. Where we differ from traditional approaches is in our definition of the field correlation time $\tau _c$. In particular, we use the definition often encountered in the literature [35–37] We emphasize that this definition of $\tau _c$ does not rely on the exact form of $g_1(\tau )$, and hence does not require numerical fitting to any particular model for $g_1(\tau )$ that must be chosen in advance, as in traditional LSCI approaches. For reference, commonly used models are $g_1(\tau ) = e^{-(\tau /\hat {\tau }_c)^n}$ for $n$ = 0.5, 1 or 2. We find that, for these commonly used models, $\tau _c = \hat {\tau }_c$ for $n$ = 0.5 or 1, and $\tau _c = \sqrt {\pi /2} \hat {\tau }_c$ for $n$ = 2.2.2 Link between correlation time and speckle contrast
In practice, when performing LSCI, not only must spatial filtering be taken into account (in the form of $\beta$) but also temporal filtering owing to the finite response time (exposure time $T$) of the camera. That is, the measured intensity at any pixel is given by $\hat {I}(t) = \int _{-\infty }^{\infty } F\left (t-t'\right ) I(t^{\prime })d t^{\prime }$, where the camera response filter is given by
A key premise of LSCI is that $\tau _c$ can be estimated from the measured speckle contrast. To this end, we define a generalized contrast based on the time-lagged covariance of the measured intensitiesCombining Eqs. (1)–5, we find that $K(\tau )$ and $g_1(\tau )$ are linked by the convolution [36]
where $Q(\tau )=\int _{-\infty }^{\infty } F\left (t^{\prime }+\tau \right ) F\left (t^{\prime }\right )d t^{\prime }$, given here byFinally, we arrive at the link between measured speckle contrast and correlation time $\tau _c$, established by
We emphasize again that this link is general and makes no recourse to a specific model for $g_{1}(\tau )$. Remarkably, this link is independent of camera exposure time (overlooking the issue of noise for now, which prescribes an optimal exposure time – see Appendix A), meaning that it can be implemented even with cameras of modest frame rate $R$. However, we note that Eq. (8) is more formal than practical. In practice, we do not have access to $K(\tau )$ over continuous time lags but instead only over discrete time lags given by the inverse frame rate $1/R$. We assume here $T \approx R^{-1}$. That is, the integral in Eq. (8) cannot be evaluated exactly but rather must be estimated from the discrete values $K_m=K(m T)$, $m=0,1,..,M$, where $M$ is assumed large enough that $K_{M}^2\xrightarrow {}0$. To perform this estimate, we use the trapezoidal rule [39]
Equation 9 represents the fundamental link between $K$ and $\tau _c$ that we will use throughout this paper. Again, this link is independent of camera exposure time. However, for this independence to remain valid even for arbitrarily small $\tau _c$ it is important that the trapezoidal rule be used to estimate the integral in Eq. (9) rather than another rule (e.g. Simpson’s rule). The trapezoidal rule makes the approximation that the $K_m$ values are pairwise connected by linear segments, meaning that, in the case $\tau _c \ll T$, the decay from $K_0$ to $K_1 \simeq 0$ becomes triangular, in accord with the shape of $Q(\tau )$. This leads to
in agreement with the simplified relation between $K^2$ and $\tau _c$ in [27]. Note that Eq. (9) is valid for $T \approx R^{-1}$. Appendix B considers the more general case where $T \leq R^{-1}$.Figure 1 illustrates our method of directly obtaining $\tau _c$ from $K^2$ for two commonly used $g_1(\tau )$ models. Specifically, $g_1(\tau ) = e^{-\tau /\tau _c}$ with $\tau _c = 0.1$ ms (Fig. 1(a)), and $g_1(\tau ) = e^{-\sqrt {\tau /\tau _c}}$ with $\tau _c = 0.5$ ms (Fig. 1(b)), represent typical dynamics in medium-sized blood vessels and parenchyma respectively [32]. We illustrate our method in the cases of two very disparate camera exposure times $T = 0.3$ ms and $T = 3$ ms. The sampled contrast values at discrete time lags (multiples of $T$) are used to estimate $\tau _c$ using Eq. (9). For both cases, $\tau _c$ is reliably recovered, with relative errors of only $0.007\%$ and $0.003\%$ for Fig. 1(a) and (b), despite the fact that the exposure times are considerably longer than $\tau _c$, even by an order of magnitude.
Fig. 1. Simulated reconstruction of $\tau _c$ based on Eq. (8) for different forms of $g_1(\tau )$. (a) $g_1(\tau ) = e^{-\tau /\tau _c}$ with $\tau _c = 0.1$ms; (b) $g_1(\tau ) = e^{-\sqrt {\tau /\tau _c}}$ with $\tau _c = 0.5$ ms ($\beta = 1$). The contrast $K^2(\tau )$ is shown sampled with different camera exposure times $T = 0.3$ ms and $T = 3$ ms.
2.3 Extension to non-ergodic measurements
So far, our approach for evaluating speckle correlation time in a model free manner has rested on the key assumption that our measurements are ergodic. That is, we made the assumption that our overall measurement time was long enough that all of statistical space was explored, meaning that temporal averaging was tantamount to ensemble averaging. We consider now regimes where this assumption breaks down. In particular, under the presence of static scattering, our measurements can no longer be considered ergodic, and the Siegert relation (Eq. (1)) can no longer be applied [19,40]. The consequences of static scattering in LSCI have been comprehensively studied [21,25,41,42]. Here, we follow conventional theoretical analyses with the aim of incorporating the possibility of non-ergodic measurements into our method of evaluating correlation times without numerical fitting.
We start with the well-known modification that must be made to the Siegert relation in the presence of static scattering. The intensity autocorrelation function $g_2(\tau )$ becomes
where $\beta$ is again the correction factor that accounts for spatial filtering only. Here, $\xi _0 = \frac {I_s}{I_d+I_s}$ and $\xi _1 = 1-\xi _0 = \frac {I_d}{I_d+I_s}$ are the fractions of dynamic and static scattering, where $I_s$ is the averaged intensity of static scattered light, and $I_d$ is the averaged intensity of dynamic scattered light with field autocorrelation function $g_1(\tau )$. A determination of $\xi _0$ is crucial here, and can be obtained in many ways. Here, we use the method prescribed in Ref. [41], which evaluates $\xi _0$ based on the discrepancy between the speckle contrast measured spatially (presumed ergodic) and measured temporally (possibly non-ergodic), using the relationOnce a value for $\xi _0$ is established, temporal contrast $K(\tau )$ is related to $g_1(\tau )$ by incorporating the modified Siegert relation (Eq. (11)), obtaining
Similar calculations as performed in Eq. (8) lead to
The final expression for $\tau _c$, allowing for the presence of static scattering, is then
3. Methods
3.1 System design and processing
Figure 2 shows our experimental setup. In brief, the illumination light is delivered by a HeNe laser (17 mW, N-LHP-925, Newport) and expanded to span the sample surface. The backscattered light is collected by an objective lens (10$\times$, 0.3NA, Olympus PLN) and 200 mm tube lens, and sent through a crossed linear polarizer before being imaged by the camera (acA2040-180km/kc, Basler).
Fig. 2. Procedure for direct reconstruction of correlation time. (a) Experimental layout. An expanded HeNe laser beam uniformly illuminates the tissue (e.g. exposed mouse brain). The backscattered speckle is collected through an objective and tube lens L, and detected by a camera (a crossed polarizer P reduces the contribution from surface scattering). A total of $N$ sequential images is recorded with exposure time $T$ and frame rate $R$. (b) For $N$ raw images, the contrast $K$ associated with different time lags $\tau _r$ is calculated, with $\tau _r$ being a multiple of the frame interval. The contrast stack is pixel-wise temporally integrated and scaled by $2/\beta$ to generate a correlation time map $\tau _c$. (c) and (d) are example intensity traces obtained from different tissue regions in b (blue square in parenchyma; red square in blood vessel). Insets show example field autocorrelation $g_1(\tau )$.
Figure 2 shows our experimental setup. In brief, the illumination light is delivered by a HeNe laser (17 mW, N-LHP-925, Newport) and expanded to span the sample surface. The backscattered light is collected by an objective lens (10$\times$, 0.3NA, Olympus PLN) and 200 mm tube lens, and sent through a crossed linear polarizer before being imaged by the camera (acA2040-180km/kc, Basler).
An imaging session consists in acquiring $N$ sequential frames recorded with camera exposure time $T$ and frame rate $R \approx T^{-1}$. The time lag $\tau _r$ is given by integer multiples of the frame interval $1/R$. In accord with the definition in Eq. (5), $K^2(\tau _r)$ is calculated from the normalized covariance of the intensity with increasing time lags (Fig. 2(b)). Then, for each pixel, the integral of $K^2(\tau _r)$ over all time lags is estimated from Eq. (9), providing a direct reconstruction of $\tau _c$.
Two important steps must be taken beforehand. First, the spatial filtering parameter $\beta$ must be pre-calibrated. This is done only once, prior to all imaging experiments, using a static calibration sample (a teflon block). In our system, $\beta = 0.49$. Note that $\beta$ here characterizes spatial filtering only, due to the camera pixel size, and is thus different from $\beta _t$ in DLSI [32] which also incorporates temporal filtering and can be spatially varying. Second, for the arguments leading to Eq. (9) to be valid, we must establish that our measurements are indeed ergodic in nature. For this we compare spatial and temporal contrast according to Eq. (12). As shown below, we consider both ergodic and non-ergodic regimes.
3.2 Phantom imaging
A liquid phantom, intended to mimic dynamic biological tissue, was made of diluted milk in a solid polydimethylsiloxane (PDMS) mold on a microscope slide. A single capillary, intended to mimic a blood vessel, of $250 \mu$m diameter was embedded in the phantom with milk pumped through. The flow speed was controlled at 1 $\mu$L/min by a syringe pump (New Era Pump Systems, Inc). For Fig. 3, $N = 1000$ images were acquired with $T = 5$ ms, $R = 197.5$ FPS to ensure that temporal average equals ensemble average.
Fig. 3. Validation experiment with a dynamic liquid phantom. (a) Schematic of the liquid phantom. (b) Stack of contrast maps $K^2(\tau _r)$ evaluated at increasing time lags $\tau _r$. (c) Decay traces of $K^2(\tau _r)$ averaged over two different ROIs, corresponding to blue (background) and red (capillary) squares in (b). (d) Inverse correlation time (ICT) map reconstructed with our direct method compared to conventional fitting methods. From top to bottom: our method using direct integration of $K^2/\beta$; method of single-exposure contrast fitting to the $g_1(\tau ) = \text {exp}(-\tau /\tau _c)$ model (MO/SU); method of single-exposure contrast fitting to the $g_1(\tau ) = \text {exp}(-\sqrt {\tau /\tau _c})$ model (MU); our method with increased camera exposure time ($T\gg \tau _c$, corresponding here to a frame rate of 17 Hz). (e) Line profiles along the dashed line in the ICT maps in (d). Scale bars in (d) and (e): $50 \mu$m.
The solid phantom, intended to mimic static biological tissue, was made of a mixture of PDMS and TiO$_2$ powder (1mg TiO$_2$ per gram of PDMS) which provides static scattering. The solid phantom was traversed by a single capillary of $100 \mu$m diameter which was connected to the syringe pump. The capillary was $220\mu$m below the scattering phantom surface. For Fig. 4, the syringe pump forced diluted milk at speeds ranging from 1$\mu$L/min to 6$\mu$L/min. At each pump speed, $N = 400$ total images were acquired with $T = 1$ ms and $R = 378$ FPS. The contrast stack $K(\tau _r)$ was calculated with $\tau _r = \frac {m}{R}, m = 0, 1,\ldots 60$ at each pixel and then averaged in a $7\times 7$ window. The spatial contrast was calculated within each $7\times 7$ pixel window and averaged across 400 frames.
Fig. 4. Validation of direct reconstruction approach for $\tau _{c}$ using a phantom with static background. (a) Milk flows through a capillary embedded in a solid scattering medium. (b) resultant ICT ($\tau _{c}^{-1}$) obtained from Eq. (9) without correction for static scattering. (c) resultant ICT ($\tau _{c}^{-1}$) obtained from Eq. (16) with correction for static scattering. (d) Corresponding ICT values obtained by direct reconstruction and a traditional fitting-based approach for different flow speeds. A failure to correct for static scattering leads to an underestimate of the ICT. Results obtained with our direct approach and a traditional fitting-based approach are consistent.
3.3 Mouse preparation and imaging
All procedures were approved by the Institutional Animal Care and Use Committee (IACUC) at Boston University, and practices were consistent with the Guide for the Care and Use of Laboratory Animals and the Animal Welfare Act. The animals used were adult male wild type mice (C57BL/6J, $\#$0000664, The Jackson Laboratory, Bar Harbor, Maine), group housed and kept on a 12-hr light/dark schedule in a vivarium with food and water available ad libitum. During surgery, anesthesia was induced and maintained using mask inhalation (1-1.5$\%$ [v/v] isoflurane–oxygen mixture, Henry Schein, Melville, New York). Animals were injected with buprenorphine (3.25 mg/kg; SC, Patterson Veterinary, Greeley, Colorado), meloxicam (5 mg/kg; SC, Covetrus, Dublin, Ohio), and dexamethasone (Henry Schein, 2.5 mg/kg; SC) and head-fixed using a stereotax (Kopf Instruments, Tujunga, California). The craniotomy was performed on a 3 mm diameter circle measured over somatosensory cortex. A sterilized 5 mm diameter glass coverslip was placed gently on top of the craniotomy and adhered to the skull with super glue (cyanoacrylate, Krazy glue, High Point, North Carolina). A thin layer of super glue was applied to any exposed skull to improve adherence of the headplate to the skull. The headplate (Scientific Instrument Facility, Boston University, http://sif.bu.edu) was then attached to the skull surrounding the window using dental cement (Ortho-Jet Liquid, Black, Lang Dental, Wheeling, Illinois).
The dataset for Fig. 5 contains 4000 total frames taken at $T = 0.8$ ms, $R = 1160$ FPS with an image size of $300\times 1000$ pixels. For both our direct reconstruction and single-exposure LSCI approaches, 1000 frames were used for processing. For the MESI approach, images with different exposures were synthesized from 4000 frames in total. Synthesized exposure times obtained by frame binning ranging from $0.8$ ms to $32$ ms were then used to produce different contrast maps based on numerical fitting. The dataset for Fig. 6 and Fig. 7 contains 3000 total frames taken at $T = 0.6$ ms, $R = 1510$ FPS with an image size of $200\times 1000$ pixels and 200 total frames taken at $T = 6$ ms, $R = 165$ FPS of image size $1000\times 1000$ pixels, respectively. For all mouse imaging experiments, our direct reconstruction was performed using contrast maps calculated with $\tau _r = \frac {m}{R}, m = 0, 1,\ldots 30$.
Fig. 5. Measurements of inverse correlation time in mouse brain. (a) ICT map processed using our method compared to conventional methods. From top to bottom: our method using direct integration of $K^2/\beta$; method of single-exposure contrast fitting to the $g_1(\tau ) = \text {exp}(-\tau /\tau _c)$ model (MO/SU); method of single-exposure contrast fitting to the $g_1(\tau ) = \text {exp}(-\sqrt {\tau /\tau _c})$ model (MU); MESI method of multiple- exposure contrast fitting to the $g_1(\tau ) = \text {exp}(-\tau /\tau _c)$ model (MO/SU). (b) Line profiles along dashed lines in the ICT maps in (a). Scale bars in (a) and (b): $50 \mu$m.
Fig. 6. Measurement of inverse correlation time with different camera exposure times and frame rates. ICT map evaluated with $T = 0.6$ ms, $R = 1510$ FPS (a), $T = 1.8$ ms, $R = 503$ FPS (b), $T = 3.6$ ms, $R = 252$ FPS (c) and $T = 4.8$ ms, $R = 47$ FPS (d). Note that $TR\approx 1$ for (a-c) and $TR\approx 0.2$ for (d). (e) Comparison of ICT values selected from random pixels in (a) and (c). The black dashed line has a slope of 1. Scale bar: 50 $\mu$m.
Fig. 7. Correction for static scattering. (a) Spatial contrast map $K_s^2$. Inset: Fraction of scattered light that is static. The circled region highlights an artifact caused by the presence of static scattering. (b) ICT ($\tau _c^{-1}$) before correction for static scattering, and (c) after correction. Scale bars in (b) and (c): $50 \mu$m.
4. Results
To verify the accuracy and robustness of our direct $\tau _c$ reconstruction approach, we imaged both calibrated phantoms (ergodic and non-ergodic) and in-vivo mouse brains, and systematically compared our results with traditional reconstruction approaches based on numerical fitting.
4.1 Phantom imaging
4.1.1 Ergodic dynamic background
An ergodic measurement was first considered where a capillary was embedded in a dynamically scattering liquid background (diluted milk), as described above (Fig. 3(a)). Directional flow in the capillary was generated by forcing milk through the capillary at a fixed speed. The total observation time was long enough ($\sim$5s) to ensure that our measurements could safely be considered ergodic. Figure 3(b) shows an example contrast stack of increasing $\tau _r$ for two regions of interest (ROIs) highlighted in blue (background) and red (capillary). The corresponding $K^2/\beta$ traces obtained from these regions are plotted in Fig. 3(c), where the blue dashed trace (background) is observed to decay smoothly over several time points, while the red dashed trace (capillary) exhibits an almost complete decay to zero within the first time point.
In Fig. 3 we compare the inverse correlation time (ICT) maps reconstructed from our method as well as from traditional methods. The top ICT map is from our direct integration method (Eq. (9)). The second and third maps are from the single-exposure LSCI method described in Ref. [41], making use of different models for $g_1(\tau )$. Specifically, the second ICT map is obtained from a fit to the MO/SU model [$g_1(\tau ) = \text {exp}(-\tau /\tau _c)$], while the third is obtained from a fit to the MU model [$g_1(\tau ) = \text {exp}(-\sqrt {\tau /\tau _c})$] [30,43]. The top three ICT maps are derived from the same dataset obtained with a moderate exposure time $T = 5$ ms. The line profile is shown in Fig. 3(e). We observe that our result matches the result of fitting to $g_1(\tau ) = \text {exp}(-\sqrt {\tau /\tau _c})$ in the background region where the scatterers exhibit diffusive motion, and matches the result of fitting to $g_1(\tau ) = \text {exp}(-\tau /\tau _c)$ in the capillary region where there is directional flow. In addition, the direct $1/\tau _c$ captures the smooth transition between two regions where mixed dynamics are present. We confirm that our direct integration methods can simultaneously capture the correct time dynamics for both regions without recourse to any model for $g_1(\tau )$ in the ergodic measurements, while more traditional methods capture the correct dynamics only for one region or another, depending on the selected model.
For comparison, we also show a fourth ICT map derived from images taken at $T = 60$ ms. In this case $T\gg \tau _c$ for all dynamics in the sample, and the contrast converges to $K^2 = \beta T/\tau _c$ [27], enabling the direct reconstruction of the ICT from only the first time point. As can be seen, our direct reconstruction procedure leads to similar ICT maps regardless of camera exposure time, even for exposure times differing by more than an order of magnitude.
4.1.2 Non-ergodic static background
We next performed validation measurements with a phantom sample in a non-ergodic regime. For this, we imaged a microfluidic capillary embedded within a static scattering background, as described above (Fig. 4(a)). Again, to mimic blood flow, milk was forced through the capillary with a syringe pump at increasing speeds.
Figure 4(b) and (c) compares the ICT maps calculated without and with correction for static scattering. In this case, if ergodicity is assumed and Eq. (9 is used to evaluate $\tau _c$, the influence of the static background leads to an underestimate of the ICT. To correct for the influence of the static background, we use Eq. (16) instead, leading to a corrected ICT map (Fig. 4(c)) and a corrected higher flow speed. Our direct reconstruction results are compared with those obtained with the traditional LSCI approach [41] using same raw data (Fig. 4(d)). In particular, this traditional approach also makes allowances for static scattering, but requires numerical fitting to a model which, because we are interested in flow dynamics, is assumed to be $g_1(\tau ) = \text {exp}(-\tau /\tau _c)$. We find that, after correction for static scattering, both our direct reconstruction approach and the traditional fitting approach yield similar ICT results. This consistency confirms the robustness of our direct reconstruction approach, which, we emphasize again, made no recourse to numerical fitting.
4.2 Mouse brain imaging
We next demonstrate our direct reconstruction method by performing vasculature imaging in the mouse brain in-vivo. For the example presented here, the brain tissue motion was determined to be sufficiently dynamic that our measurements could be considered roughly ergodic (the static contribution $\xi _0$ was determined to be less than 5% on average, and the error that comes from ignoring 5% static scattering with $\eta _1 = \sqrt {2}$ is only about 3.5% according to Eq. (16)). Figure 5 shows a comparison of the resultant ICT values obtained from our direct method and from traditional model-based fitting methods. In Fig. 5(a), the top three ICT maps are from direct reconstruction and from single exposure LSCI with numerical fitting to the MO/SU and MU models respectively. The fourth (bottom) ICT map is obtained from the multiple exposure speckle imaging (MESI) method described in Ref. [21], making use of a fit to the MO/SU model and therefore similar to the second single-exposure result. The corresponding line profiles from these ICT maps are shown in Fig. 5(b), in which the vessel and parenchyma regions are highlighted by arrows. We observe that our direct reconstruction method simultaneously matches the results from the MO/SU model when associated with blood-vessel regions, and the results from the MU model when associated with parenchymal regions [27,34], demonstrating that our method is able to capture the correct dynamics in either region of interest, without any assumption of a model or any requirement of model fitting (in contrast to traditional methods).
We further verify that our direct reconstruction method provides consistent results regardless of the camera exposure time. Figure 6 shows an example where a small field of view was imaged at $T = 0.6$ ms, 1510 FPS. The raw frames were binned by factors of 1, 3 and 6 to synthesize increasing exposure durations ($T = 0.6, 1.8$ and 3.6 ms and $T\approx 1/R$). The corresponding ICT maps are shown in Fig. 6(a-c). We arbitrarily selected 1000 pixels and plotted the ICT value for each pixel evaluated at $T = 0.6$ ms and $T = 3.6$ ms as shown in Fig. 6(e). The scatter plot reveals a good match between the two ICT maps, confirming the relative insensitivity of our direct reconstruction method to exposure time. In addition, we analyzed images from raw frames binned by a factor of 8 and down-sampled by a factor of 4, resulting in an effective exposure time $T = 4.8$ ms and frame rate $R = 47$ FPS (i.e. $T \neq 1/R$). In this case, one has to account for the discrepancy between exposure time and frame period in our evaluation of $K_0^2$ (Appendix B). The corresponding ICT map calculated from Eq. (21) is shown in Fig. 6(d), where we continue to capture the tissue dynamics even at this much slower near-video frame rate. We note, however, that in practice considerations of noise and statistical sampling prescribe an optimal exposure time and frame rate, as discussed in Appendices A and B.
Finally, we demonstrate the application of our method to non-ergodic in-vivo measurements, which often occur when performing mouse-brain imaging [27,29]. Figure 7(a) illustrates an example where we observed the presence of a localized static field, apparent as a bright artifact (circled region) in the spatial contrast image, and also apparent in the static fraction map $\xi _0$ (Fig. 7(a) inset). The resultant ICT ($\tau _c^{-1}$) maps obtained without (Eq. (9)) and with (Eq. (16)) static-field correction are shown in Figs. 7(b) and 7(c) respectively. We observe that without correction, the resultant ICT is underestimated in the presence of static scattering.
5. Discussion
In summary, we present a temporal LSCI approach to extract correlation times from measurements of speckle contrast that differs from conventional approaches in several respects. To begin, we make use of a generalized definition of correlation time that, in the case of ergodic measurements, altogether circumvents the problem of having to predefine a theoretical model for $g_1(\tau )$ governing the speckle dynamics. This has several advantages. Because our reconstruction method is model-agnostic, it does not require adjusting a model to different image regions, as in [32]. Moreover, our reconstruction method does not require computationally intensive numerical fitting to any particular model but is performed instead by simple summation of contrast images (Eq. (9)), with processing times thus greatly reduced (see Appendix C). Finally, our method remains valid even if the camera is only of modest speed, meaning that a specialized high-speed camera is not required and larger fields of view can be addressed. For comparison, an alternative Fourier transform laser speckle imaging (FT-LSI) method has been proposed that also provides quantitative reconstructions based on Fast Fourier Transforms [43]. However, this method requires a fast camera ($\sim$ 10KHz) to capture the full speckle-fluctuation power spectral density. We note that while in theory our approach is independent of camera speed, in practice an optimized camera speed is nevertheless prescribed when noise is taken in to account. In general, our approach remains very tolerant to wide variations in speed, and, indeed, we have demonstrated that it remains accurate when the speckle dynamics fully decay within even a single exposure time.
For added versatility we additionally provide an extension of our approach to allow for the possibility of measurements that are not fully ergodic. In other words, we allow for the possibility that a portion of the measured speckle field is static. Our method loses some of its attractiveness in this case in that it is no longer strictly model free. However, only the general shape of $g_1(\tau )$ is required in this case, in the form of a single parameter $\eta _1$, rather than any detailed model for $g_1(\tau )$. Moreover, our method continues to benefit from its ability to reconstruct correlation times without recourse to any numerical fitting, thus maintaining reduced processing times.
In a practical experiment the temporal resolution has to be taken into consideration. For $M$ time-lagged contrast measurements acquired over a total observation time $T_M = N_0MT$, $MT$ should be chosen long enough such that $K^2(MT)$ fully decays to zero. Moreover, the total observation time $T_M$ should be chosen substantially longer than $MT$ to ensure that the intensity covariance at different time lags explores the full statistics. In our case $MT$ was at most 20 ms and we recorded over an observation time $T_M = 1\sim 2$ s, such that $N_0$ was greater than 100, which effectively limited the statistical error to less than 10$\%$ [16]. Note that this total observation time is similar to that used with DLSI [32], but without the requirement of ultrafast frame rates.
In closing, we note that our direct correlation-time reconstruction approach can be applied to other measurement modalities involving speckle dynamics, such as diffusive wave spectroscopy [44], speckle visibility spectroscopy [13] or dynamic light scattering imaging [27,32]. The simplicity of our approach can also be attractive for clinical applications, including the monitoring of cerebral blood flow [7,45], retinal flow [46,47] and skin perfusion [4,48], to name a few.
Appendix A: Optimum camera exposure time
When evaluating the reliability of our approach to calculate $\tau _c$, we must consider noise. Because $K^2(\tau )$ decreases with increasing $\tau$, its integral (Eq. (9)) is most influenced by the first few values of $K^2_m$. We thus focus on evaluating the signal to noise ratio (SNR) associated with $K^2_0=K^2(0)$ (indeed, when $\tau _c \ll T$, the integral is derived entirely from $K^2_0$). Several noise sources affect the measure of $K^2_0$. For example, shot-noise and camera readout noise produce a bias in $K^2_0$, which should be corrected for. As shown in [36], the bias is given by
However, the above analysis does not take into account uncertainties due to insufficient sampling, which are often the dominant source of noise. In general, the discrepancy between the unbiased sample variance and the true variance $\sigma ^2$ of a random variable is given by [49]
where $N$ is the number of samples. This discrepancy (also discussed in [20]) leads to an error in the estimated contrast given by (for $N \gg 1$)For a total measurement time $T_M$, we have $N \approx T_M/T$. We thus find that the relative error scales as $N^{-1} \approx T/T_M$, suggesting that shorter exposure times are better.
Taking into account both detection noises and sampling error, we can very roughly derive an optimum $T$ by maximizing the ratio $K^2_0/(K^2_{bias} + \delta K^2_0)$. This leads to
For example, using parameters from our experiments $\sigma _n = 12$ electrons, $\left\langle I \right\rangle \approx 3.2\times 10^5$ photoelectrons/s, $T_M = 1$ s, $\beta = 0.49$ and $\tau _c \approx 0.2$ ms, we obtain $T_{opt}\sim 4$ ms. One must be careful with this result since we have treated $K^2_{bias}$ as a noise term rather than the uncertainty in $K^2_{bias}$. Nevertheless, Eq. (20) indicates that $T$, in the presence of inevitable readout noise, should be neither too short nor too long. Note that our result for $T_{opt}$ is in accord with the phenomenologically observed optimal exposure time reported in [50].
Appendix B: Case when $T\leq 1/R$
In the situation where the camera exposure time $T$ is smaller than the frame duration $1/R$ (i.e. the filtering and sampling times are different), one needs to correct for the discrepancy introduced in the zero-time lagged contrast; that is, $K_0$ is only affected by the integration time $T$ whereas $K_m^2, m = 1,2,{\ldots }M$ is sampled at each $\tau _r = m/R$. In this case, we use the geometry shown in Fig. 8 to estimate the integration of the first interval separately, leading to
Fig. 8. Estimate of the integration of $K^2(\tau )$ with discrete data points $K_m^2, m = 0,1,{\ldots }M$. The first interval is sub-divided into two parts (a triangle and a trapezoid). The height of the trapezoid at time zero is estimated by interpolation using $K_1^2$ and $K_2^2$. The area of each component is then estimated geometrically and summed.
We note that Eq. (21) reduces to Eq. (9) when $T = 1/R$, and indeed, in this case both equations provide near exact measures of $\tau _c$. In the event $T < 1/R$, Eq. (21) generally provides a more accurate estimate of $\tau _c$ than Eq. (9), but with an error that increases with increasing discrepancy between $T$ and $1/R$, though remaining small if $T$ is much greater than $\tau _c$. For example, $TR=0.2$ for Fig. 6(d), and the relative error in our estimate of $\tau _c$ was less than 5% for the blood vessel and 20% for the parenchyma. In general, it is preferable to operate with $TR$ close to 1.
Appendix C: Comparison of processing times
Table 1 shows typical processing times for a $300\times 1000 \times 400$ image stack using three methods. The processing is performed on a desktop computer equipped with a i7-6700K CPU and 64GB RAM. For the direct reconstruction method used here, 30 contrast maps $K^2(\tau )$ are calculated from a total of 400 raw intensity images. The correlation time is derived from the direct integration of $K^2(\tau )$ at each pixel. For the method of single exposure LSCI, a spatial contrast map with a $7\times 7$ moving window size and a temporal contrast map with $N_t = 25$ temporal window size are calculated. The two contrast maps are used to find the fraction of static scattering and then fitted to a theoretical model for speckle decorrelation following procedures in [41]. For the method of MESI, a total of 3000 raw images are used. The images are binned to synthesize images at 30 different exposure times, from which 30 contrast maps are calculated for numerical fitting [21].
Table 1. Processing time comparison of the direct reconstruction method versus single-exposure LSCI and MESI
Funding
National Institutes of Health (R01EB029171).
Acknowledgments
We thank Lisa Kretsge, Kelly Wingfield, and Alberto Cruz-Martin for providing mouse samples. We thank David Boas for helpful discussions and suggestions.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. A. Fercher and J. D. Briers, “Flow visualization by means of single-exposure speckle photography,” Opt. Commun. 37(5), 326–330 (1981). [CrossRef]
2. J. Senarathna, A. Rege, N. Li, and N. V. Thakor, “Laser speckle contrast imaging: theory, instrumentation and applications,” IEEE Rev. Biomed. Eng. 6, 99–110 (2013). [CrossRef]
3. P. G. Vaz, A. Humeau-Heurtier, E. Figueiras, C. Correia, and J. Cardoso, “Laser speckle imaging to monitor microvascular blood flow: a review,” IEEE Rev. Biomed. Eng. 9, 106–120 (2016). [CrossRef]
4. D. A. Boas and A. K. Dunn, “Laser speckle contrast imaging in biomedical optics,” J. Biomed. Opt. 15(1), 011109 (2010). [CrossRef]
5. D. Briers, D. D. Duncan, E. R. Hirst, S. J. Kirkpatrick, M. Larsson, W. Steenbergen, T. Stromberg, and O. B. Thompson, “Laser speckle contrast imaging: theoretical and practical limitations,” J. Biomed. Opt. 18(6), 066018 (2013). [CrossRef]
6. D. A. Boas, Diffuse Photon Probes of Structural and Dynamical Properties of Turbid Media: Theory and Biomedical Applications (University of Pennsylvania, 1996).
7. A. K. Dunn, H. Bolay, M. A. Moskowitz, and D. A. Boas, “Dynamic imaging of cerebral blood flow using laser speckle,” J. Cereb. Blood Flow Metab. 21(3), 195–201 (2001). [CrossRef]
8. H. Cheng, Q. Luo, S. Zeng, S. Chen, J. Cen, and H. Gong, “Modified laser speckle imaging method with improved spatial resolution,” J. Biomed. Opt. 8(3), 559–564 (2003). [CrossRef]
9. P. Li, S. Ni, L. Zhang, S. Zeng, and Q. Luo, “Imaging cerebral blood flow through the intact rat skull with temporal laser speckle imaging,” Opt. Lett. 31(12), 1824–1826 (2006). [CrossRef]
10. O. Thompson, M. Andrews, and E. Hirst, “Correction for spatial averaging in laser speckle contrast analysis,” Biomed. Opt. Express 2(4), 1021–1029 (2011). [CrossRef]
11. J. Qiu, Y. Li, Q. Huang, Y. Wang, and P. Li, “Correcting speckle contrast at small speckle size to enhance signal to noise ratio for laser speckle contrast imaging,” Opt. Express 21(23), 28902–28913 (2013). [CrossRef]
12. J. Ramirez-San-Juan, R. Ramos-Garcia, G. Martinez-Niconoff, and B. Choi, “Simple correction factor for laser speckle imaging of flow dynamics,” Opt. Lett. 39(3), 678–681 (2014). [CrossRef]
13. R. Bandyopadhyay, A. Gittings, S. Suh, P. Dixon, and D. J. Durian, “Speckle-visibility spectroscopy: A tool to study time-varying dynamics,” Rev. Sci. Instrum. 76(9), 093110 (2005). [CrossRef]
14. C. P. Valdes, H. M. Varma, A. K. Kristoffersen, T. Dragojevic, J. P. Culver, and T. Durduran, “Speckle contrast optical spectroscopy, a non-invasive, diffuse optical method for measuring microvascular blood flow in tissue,” Biomed. Opt. Express 5(8), 2769–2784 (2014). [CrossRef]
15. C. Huang, D. Irwin, Y. Lin, Y. Shang, L. He, W. Kong, J. Luo, and G. Yu, “Speckle contrast diffuse correlation tomography of complex turbid medium flow,” Med. Phys. 42(7), 4000–4006 (2015). [CrossRef]
16. A. C. Völker, P. Zakharov, B. Weber, F. Buck, and F. Scheffold, “Laser speckle imaging with an active noise reduction scheme,” Opt. Express 13(24), 9782–9787 (2005). [CrossRef]
17. D. D. Duncan, S. J. Kirkpatrick, and R. K. Wang, “Statistics of local speckle contrast,” J. Opt. Soc. Am. A 25(1), 9–15 (2008). [CrossRef]
18. S. J. Kirkpatrick, D. D. Duncan, and E. M. Wells-Gray, “Detrimental effects of speckle-pixel size matching in laser speckle contrast imaging,” Opt. Lett. 33(24), 2886–2888 (2008). [CrossRef]
19. J.-Z. Xue, D. Pine, S. T. Milner, X.-L. Wu, and P. Chaikin, “Nonergodicity and light scattering from polymer gels,” Phys. Rev. A 46(10), 6550–6563 (1992). [CrossRef]
20. P. Zakharov, “Ergodic and non-ergodic regimes in temporal laser speckle imaging,” Opt. Lett. 42(12), 2299–2301 (2017). [CrossRef]
21. A. B. Parthasarathy, W. J. Tom, A. Gopal, X. Zhang, and A. K. Dunn, “Robust flow measurement with multi-exposure speckle imaging,” Opt. Express 16(3), 1975–1989 (2008). [CrossRef]
22. P.-A. Lemieux and D. Durian, “Investigating non-gaussian scattering processes by using nth-order intensity correlation functions,” J. Opt. Soc. Am. A 16(7), 1651–1664 (1999). [CrossRef]
23. D. A. Boas and A. G. Yodh, “Spatially varying dynamical properties of turbid media probed with diffusing temporal light correlation,” J. Opt. Soc. Am. A 14(1), 192–215 (1997). [CrossRef]
24. C. Ayata, A. K. Dunn, Y. Gursoy-Özdemir, Z. Huang, D. A. Boas, and M. A. Moskowitz, “Laser speckle flowmetry for the study of cerebrovascular physiology in normal and ischemic mouse cortex,” J. Cereb. Blood Flow Metab. 24(7), 744–755 (2004). [CrossRef]
25. P. Zakharov, A. Völker, M. Wyss, F. Haiss, N. Calcinaghi, C. Zunzunegui, A. Buck, F. Scheffold, and B. Weber, “Dynamic laser speckle imaging of cerebral blood flow,” Opt. Express 17(16), 13904–13917 (2009). [CrossRef]
26. J. Xu, A. K. Jahromi, J. Brake, J. E. Robinson, and C. Yang, “Interferometric speckle visibility spectroscopy (isvs) for human cerebral blood flow monitoring,” APL Photonics 5(12), 126102 (2020). [CrossRef]
27. C. Liu, K. Kılıç, S. E. Erdener, D. A. Boas, and D. D. Postnov, “Choosing a model for laser speckle contrast imaging,” Biomed. Opt. Express 12(6), 3571–3583 (2021). [CrossRef]
28. D. D. Duncan and S. J. Kirkpatrick, “Can laser speckle flowmetry be made a quantitative tool?” J. Opt. Soc. Am. A 25(8), 2088–2094 (2008). [CrossRef]
29. A. B. Parthasarathy, S. S. Kazmi, and A. K. Dunn, “Quantitative imaging of ischemic stroke through thinned skull in mice with multi exposure speckle imaging,” Biomed. Opt. Express 1(1), 246–259 (2010). [CrossRef]
30. D. J. Pine, D. A. Weitz, P. M. Chaikin, and E. Herbolzheimer, “Diffusing wave spectroscopy,” Phys. Rev. Lett. 60(12), 1134–1137 (1988). [CrossRef]
31. S. S. Kazmi, L. M. Richards, C. J. Schrandt, M. A. Davis, and A. K. Dunn, “Expanding applications, accuracy, and interpretation of laser speckle contrast imaging of cerebral blood flow,” J. Cereb. Blood Flow Metab. 35(7), 1076–1084 (2015). [CrossRef]
32. D. D. Postnov, J. Tang, S. E. Erdener, K. Kılıç, and D. A. Boas, “Dynamic light scattering imaging,” Sci. Adv. 6(45), eabc4628 (2020). [CrossRef]
33. M. A. Davis and A. K. Dunn, “Dynamic light scattering monte carlo: a method for simulating time-varying dynamics for ordered motion in heterogeneous media,” Opt. Express 23(13), 17145–17155 (2015). [CrossRef]
34. M. A. Davis, L. Gagnon, D. A. Boas, and A. K. Dunn, “Sensitivity of laser speckle contrast imaging to flow perturbations in the cortex,” Biomed. Opt. Express 7(3), 759–775 (2016). [CrossRef]
35. J. W. Goodman, Speckle Phenomena in Optics: Theory and Applications (Roberts and Company Publishers, 2007).
36. J. Mertz, Introduction to Optical Microscopy (Cambridge University Press, 2019).
37. R. Loudon, The Quantum Theory of Light (OUP Oxford, 2000).
38. B. J. Berne and R. Pecora, Dynamic Light Scattering: With Applications to Chemistry, Biology, and Physics (Courier Corporation, 2000).
39. J. Stoer and R. Bulirsch, Introduction to Numerical Analysis, vol. 12 (Springer Science & Business Media, 2013).
40. P. N. Pusey and W. Van Megen, “Dynamic light scattering by non-ergodic media,” Phys. A 157(2), 705–741 (1989). [CrossRef]
41. Y. Wang, D. Wen, X. Chen, Q. Huang, M. Chen, J. Lu, and P. Li, “Improving the estimation of flow speed for laser speckle imaging with single exposure time,” Opt. Lett. 42(1), 57–60 (2017). [CrossRef]
42. J. C. Ramirez-San-Juan, C. Regan, B. Coyotl-Ocelotl, and B. Choi, “Spatial versus temporal laser speckle contrast analyses in the presence of static optical scatterers,” J. Biomed. Opt. 19(10), 106009 (2014). [CrossRef]
43. J. Buijs, J. Gucht, and J. Sprakel, “Fourier transforms for fast and quantitative laser speckle imaging,” Sci. Rep. 9(1), 13279 (2019). [CrossRef]
44. J. Xu, A. K. Jahromi, and C. Yang, “Diffusing wave spectroscopy: a unified treatment on temporal sampling and speckle ensemble methods,” APL Photonics 6(1), 016105 (2021). [CrossRef]
45. A. K. Dunn, “Laser speckle contrast imaging of cerebral blood flow,” Ann. Biomed. Eng. 40(2), 367–377 (2012). [CrossRef]
46. A. Ponticorvo, D. Cardenas, A. K. Dunn, D. Ts’o, and T. Q. Duong, “Laser speckle contrast imaging of blood flow in rat retinas using an endoscope,” J. Biomed. Opt. 18(9), 090501 (2013). [CrossRef]
47. D. D. Patel and D. M. Lipinski, “Validating a low-cost laser speckle contrast imaging system as a quantitative tool for assessing retinal vascular function,” Sci. Rep. 10(1), 7177 (2020). [CrossRef]
48. B. Choi, N. M. Kang, and J. S. Nelson, “Laser speckle imaging for monitoring blood flow dynamics in the in vivo rodent dorsal skin fold model,” Microvasc. Res. 68(2), 143–146 (2004). [CrossRef]
49. S. A. Teukolsky, B. P. Flannery, W. Press, and W. Vetterling, “Numerical recipes in c,” SMR 693, 59–70 (1992).
50. S. Yuan, A. Devor, D. A. Boas, and A. K. Dunn, “Determination of optimal exposure time for imaging of blood flow changes with laser speckle contrast imaging,” Appl. Opt. 44(10), 1823–1830 (2005). [CrossRef]
