1. Introduction

Laser speckle contrast imaging (LSCI) is a relatively inexpensive widefield imaging method that allows for assessment of tissue perfusion and vascular flow across entire fields of view with high spatiotemporal resolution [1,2]. LSCI gives a measure of flow based on speckle contrast, K, or the degree of blurring of a laser speckle pattern over a spatial region in a time-integrated image due to moving scatterers. Because of the simplicity of the technique, requiring only a laser, diffuser or beam expander, optical filter, and camera, LSCI has been used extensively in neurophysiology for assessing perfusion and functional activity [3]. It has also gained some traction in clinical applications, such as in burn assessment, dermatology, and surgery [4], where it provides similar clinical capabilities to laser doppler flowmetry but with additional spatial information. LSCI, however, is limited to producing qualitative measurements of flow due to influences from the laser coherence, mismatches between the sensor pixel size and laser speckle size, and tissue scattering properties, which was addressed by the introduction of multi-exposure speckle imaging (MESI) [5–9]. In MESI, images with different exposure times are acquired and fit to a more robust model that accounts for laser coherence and mismatch between speckle size and pixel sensor size with the normalization parameter, β, for tissue scattering properties with the fraction of dynamic scattering parameter, ρ, and for directional vs. random scatterer motion by varying the model of light scattering in the tissue.

Some work has been done to bring quantitative MESI into clinical practice [10], and while the data acquisition is feasible in a clinical setting, the data analysis requires too much time to enable the real-time viewing needed for clinical application. The bottleneck for the use of whole field-of-view MESI in real-time applications is the need to fit speckle contrast values from multiple exposures at each pixel to models with the four parameters described above to reconstruct a single MESI image [11], a process that can take several minutes and up several hours depending on the image size and the machine used. As a result, datasets are typically processed by restricting analysis to smaller, sparse regions of interest (ROI) over vessels or parenchymal tissue [12]. We find that processing a 1-megapixel image takes ∼20 minutes using 24 physical cores for parallel processing. To fit every image set of that size taken over the course of a minute at a 3-Hz acquisition rate would require 60 hours, making the use of ROIs for temporal analysis a necessity. Some methods have used calibrations before or during an imaging session to estimate the β and ρ terms to improve LSCI estimations [13] or speed up MESI fitting [10,12,14], but have the limitations of not accounting for spatial or temporal variability in β, caused by differences in surface roughness [15–18], laser spectrum fluctuations [15,16], optical aberrations [17,18], or out-of-focus surfaces [17,18], or for temporal variability in ρ from, for example, changes in the fraction of tissue with flowing blood cells due to a stroke [11]. Other groups have incorporated machine learning to decrease the time required for processing MESI images [19,20], which offers vast improvements in speed but can have limitations in reliability and accuracy [21].

Here we propose a quasi-analytic solution, noted as a Rapid Estimation of Multi-exposure Imaging (REMI), with a modification to incorporate spatial averaging of some parameters denoted as spatial-REMI (sREMI), to determine the scattering model, estimate the normalization and fraction of dynamic scattering parameters, and approximate the value of the correlation time with relatively high accuracy and up to 1500 × faster than traditional least-squares fitting. The increase in speed allows for real-time assessments of perfusion for experimental or surgical applications after taking a few image sets and can be implemented with a rolling buffer of image data while omitting the need for prior static calibration. We validate the accuracy of the method using simulated data, with and without artificial noise, and compare the performance to that of least-squares fitting for simple and mixed scattering models, demonstrate further improvements in accuracy and processing time by utilization of logarithmically-spaced exposure times over the typically-used exposure times, and finally compare the in vivo measurements taken from LSCI, least-squares fitting, and sREMI in a mouse model of photothrombotic stroke.

2. Theory of laser speckle contrast imaging

For LSCI, coherent light that has reflected off a rough, randomly scattering surface is imaged on a camera sensor, where the photons interfere constructively and destructively, producing the bright-dark spots of a speckle pattern. Moving scatterers, namely cells moving through blood vessels in tissue, introduce temporal fluctuations in this interference pattern which blur the speckles within the image when integrated over a set exposure time. The degree of blurring increases with increasing speed of the moving scatterers. Speckle contrast, K, is quantified as

 

Fig. 1. Laser speckle imaging and fitting examples. (A) Raw laser speckle image (5-ms exposure, left), and processed speckle contrast image (right). (B) Processed laser speckle images at varying exposure times from 20 µs up to 500 ms. Scale bar = 200 µm. (C) Different models are fitted to the speckle contrast data taken from indicated points in A, with the red point situated on top of a larger surface vessel and the blue point situated over a parenchymal region. The top row uses the exposure times commonly found in the literature [7,32] from 50 µs – 80 ms, while the bottom row uses 20 logarithmically-spaced exposure times from 20 µs – 650 ms, with all exposures listed in Table S1. Simple-model fitting (sm-MESI) shows multiple regions of overshoots and undershoots. Mixed-model fitting (mm-MESI) and Rapid Estimation of Multi-exposure images (REMI) follow the data more closely. (D) Example MESI images produced using sm-MESI (top left), mm-MESI (top middle), REMI (bottom left), and spatial-REMI (bottom, middle). Images are 788 × 788 pixels. Required processing times are indicated below the images. 10-pixel-wide line profiles for the ROIs indicated with yellow lines are displayed for each method to the right. Background was taken as the bottom 25% of inverse correlation time values and the line profiles were normalized to this background.

Download Full Size | PDF

In practice, the electric field cannot be easily measured to characterize g₁(τ). The intensity, I, of the scattered light, however, is related to the electric field by $I = E{E^\ast }$ [1,9] and can be measured by a camera sensor. The intensity autocorrelation function, g₂(τ),

In LSCI, the β term is often determined by calibration to a static phantom. The use of a single exposure time, T, usually between 1 and 10 ms [26], allows for high temporal resolution imaging of relative blood flow speed. If the exposure is long enough compared to the anticipated correlation time, Eq. (7) can be simplified to estimate flow speed as

LSCI, however, fails to account for scattering by static components, which causes significant loss in sensitivity to moving scatterers [9,13,27]. To address the limitations of LSCI, the autocorrelation function, g₂(τ) was expanded by Parthasarathy et al [9] to include a scattering term from non-moving scatterers as

In order to fit the additional parameters, images are taken at multiple exposure times (Fig. 1(B)) with contrast values fit to Eq. (10) as a more quantitative measure in a method known as multi-exposure speckle imaging (MESI) [5–9]. The nonergodic noise normally added as a constant can be considered negligible in systems with low readout and sensor noise, or with additional noise reduction schemes [28].

Work done by Postnov, et al. [29] showed that Eq. (10) worked well for the case of multiple scattered light from scatterers undergoing ordered movement (which has the same lineshape as single scattering from unordered movement), which is seen in larger blood vessels. However, the model was less sensitive to flow changes in parenchymal tissue, where multiple scattering from scatterers undergoing unordered movement dominates [30]. Postnov, et al. thus modified g₁(τ) to account for different scattering models [25,30] with single or multiple scattering events, and for ordered or unordered movement of the scatterers. The electric field autocorrelation function was adjusted to be

3. Principles for rapid estimation of multi-exposure imaging (REMI)

Typically, regression methods seek to iteratively solve many variables simultaneously which requires many, often repeated, calculations. Our proposed method, REMI, aims to solve for the unknown variables, those being β, ρ, τ_c, and D_MU, in a specific order using a set of limiting conditions from the model proposed in Eq. (13), drastically reducing the required number of calculations. We find REMI leads to fairly optimal fits and output values nearly matching more computationally intensive least-squares methods (Fig. 1(C)) with sREMI, an extension incorporating spatial averaging of some parameters, producing images of similar clarity with fewer unrealistic spikes within the vessel profiles compared to least-squares fitting methods (Fig. 1(D)).

To begin, we consider three limiting conditions in Eq. (13): when the exposure time matches the correlation time ($T = {\tau _c}$), when the exposure time is much shorter than the correlation time ($T < < {\tau _c}$) or approaches 0, and when the exposure time is much greater than the correlation time ($T > > {\tau _c}$) or approaches ∞. For the condition where $T = {\tau _c}$, x is set to 1 in Eq. (10) and Eq. (12), which, combined in Eq. (13), gives a value, $K_{crit}^2$, as:

 

Fig. 2. Underlying principles for implementing REMI. (A) Example plot of the ordered multiple scattering (red) and unordered multiple scattering (blue) models (left) and their respective derivatives (right) with example simulated data as white squares. (B) Ratio of (K²(0)-K²(T)) to the derivative, dK²(T)/d ln(T) (solid line) for extrapolation to determine the limiting condition setting β, with a simple lowpass fit (dotted line) and the expanded fit (dashed line). (C) Ratio of (K²(T) – K²(∞)) to the derivative, dK²(T)/d ln(T) (solid line) for extrapolation to determine the limiting condition setting ρ, with a simple highpass fit (dotted line) and the expanded fit (dashed line). Solid lines are color-coded to match the scattering models from A. (D) Flowchart for the REMI algorithm. Numbered steps in the flowchart are presented in numbered bubbles in A. The gray arrows indicate the optional spatial filtering of β for sREMI.

Download Full Size | PDF

When plotting the derivatives, we see noticeable differences in peak amplitude and location of the peak in the first-order derivative with respect to correlation time for each scattering model (Fig. S2A). While the peaks appear proximal to where T = τ_c, the location for the minimum peak of dK²/d ln(T), noted as T_peak, is laterally displaced in relation to the true τ_c between ∼0.8τ_c and ∼3.3τ_c depending on the degree of dynamic scattering, ρ, and the scattering type, D_MU (Fig. S2A). Simply finding the peak of the differential will not give the true τ_c value due to the unknown displacement induced by the unresolved parameters. To solve for the missing parameters, β and ρ, we examine the difference between K²(T) and the limit as $T \to 0$ and $T \to \infty $, whereby:

4. Implementation of REMI

From Eq. (14), approximation of $K_{crit}^2$ where T = τ_c requires estimation of the scattering model, D_MU, the fraction of dynamic scattering, ρ, and the normalization factor, β. From Eq. (20) and Eq. (21), β and ρ can be estimated only when D_MU and T_peak are known. Of these, T_peak is the first variable that can be determined simply by locating the minimum peak location of dK²/d ln(T) (Fig. 1(B)). We use T_peak in Eq. (22) and Eq. (23) to approximate β and ρ assuming D_MU = 0, then compare the minimum amplitude of dK²/d ln(T) to the theoretical minimum of Eq. (16) and Eq. (17) using the estimated terms. D_MU is then set to a value between 0 and 1, determined by the distance between dK²(T_peak)/d ln(T) and the theoretical amplitudes for the two scattering models (Fig. 2(A) right, curly bracket). Once the scattering model is determined, the limiting conditions of $\mathop {\lim }\limits_{x \to 0} {K^2}$ and $\mathop {\lim }\limits_{T \to \infty } {K^2}$ are estimated by using Eq. (22) and Eq. (23), now with the appropriate D_MU, and β is set to $\mathop {\lim }\limits_{\textrm{x} \to 0} {\textrm{K}^2}$, while ρ is calculated using Eq. (15). At this point, the scattering model, D_MU, the fraction of dynamic scattering, ρ, and the normalization constant, β, have been solved. These parameters enable the use of Eq. (14) to set the $K_{crit}^2$ value which equals ${K^2}(T )$ when T = τ_c,, allowing for simple interpolation to solve. We provide more details on some of the critical steps in this process in the rest of this section.

4.1 Pre-processing discrete data

Laser speckle image data can be noisy, especially with fewer averaged speckle contrast images per exposure. Differentiating noisy discrete data only further exacerbates the noise. To reduce noise, the discrete data needs to be filtered in some way, but filtering can often suppress important higher frequency features. To balance this, we use a 1:4:1 weighted moving average filter for evenly-spaced exposures in the logarithmic domain, or a weighted, nonuniform moving average for unequally-spaced typical exposure times [32] (see Table S1; see Supplement 1 for details on nonuniform moving average filter). The central difference of the filtered data is then taken over intervals of ln(T) to achieve a semilogarithmic derivative. To further reduce noise, the derivative is filtered again using a simple 3-point moving average filter for equally spaced exposure times or a weighted-moving average filter for unequally spaced data with a width of 1 decade. Unequally spaced data is resampled using the filtered differential to interpolate and align the exposure times where the central difference is calculated for further operations.

4.2 Determining the scattering model and the unsolved parameters

To determine the scattering model, we first need to find the steepest point of the slope of the speckle contrast, K², to determine T_peak. The exposure time values of the most negative value of dK²/d ln(T) and the closest two surrounding points are used in a quadratic Newton interpolation polynomial while setting the derivative of the polynomial to 0 to locate T_peak. dK²(T_peak)/d ln(T) is determined by setting T_peak as the interpolation point of the polynomial. As a moving average filter is expected to attenuate the peak of the differential, we increase the amplitude by 5%, determined empirically, to compensate in the next 2 steps. We set γ as the ratio of T_peak to the shortest valid exposure time, T_short, for Eq. (22) and the ratio of T_peak to the longest valid exposure time, T_long, for Eq. (23). (Note, because we take the central difference of the discrete speckle contrast data, we cannot calculate dK²/d ln(T) at the absolute shortest and longest exposure times taken during acquisition, and instead use exposure times that are one step longer or shorter than the shortest or longest exposure time taken, respectively.) While assuming D_MU = 0 in Eq. Equation (20) and Eq. (21) we approximate ${K^2}(0 )$ and ${K^2}(\infty )$, respectively. Equation (1)5 is then used to solve for ρ.

To determine the correct scattering model, we compare dK²(T_peak)/d ln(T) to the theoretical max negative slope for both D_MU conditions. Repeatedly solving for these minimum values in Eq. (16) and Eq. (17) for comparison is computationally expensive. However, the peak value, when holding D_MU constant, only depends on ρ and is scaled by β. This simplification allows for empirically computing the peak values for both D_MU = 0 and D_MU = 1 as a function of ρ from 0 to 1 and fitting a polynomial for reference (this paper used a 10^th order polynomial; see Supplement 1 for polynomial expression). In this way, we can now pass an M × N matrix of ρ values and output two M × N matrices of peak values for comparison with low computational overhead. We predict β and ρ using Eq. (16) and Eq. (17) while letting D_MU = 0 and calculate the theoretical peak amplitudes for D_MU = 0 and D_MU = 1 using ρ and β with the polynomial fits. The scattering model D_MU is then given as the distance of dK²(T_peak)/d ln(T) from the theoretical peak in the D_MU = 0 condition divided by the distance between the theoretical peaks for the D_MU = 0 and D_MU = 1 conditions, bounded between 0 and 1 (Fig. 2(A); right, bubble 4).

4.3 Spatially averaging β for spatial-REMI (sREMI)

Up to this point, all data is processed in a pixel-by-pixel manner, incorporating no spatial information, which can result in significant salt-and-pepper noise (Fig. 1(D), bottom-left). The normalization constant, β, is expected to vary slowly across the field-of-view due to optical aberrations and defocus across the imaging field. Therefore, β can be spatially filtered to reduce salt-and-pepper artifacts (Fig. 1(D), bottom-right). Equation (2)2 is used with the first estimate of D_MU to calculate an initial β. A gaussian filter with a width of at least 10% of the image field is applied to the estimated β. We then re-estimate the scattering model, D_MU, as outlined above, using the spatially filtered ${K^2}(0 )$ when initially setting D_MU = 0 (Fig. 2(D), dashed outline).

4.4 Interpolation of the correlation time, τ_c

Once the scattering model has been determined, we re-evaluate ρ using Eq. (23), setting D_MU at the inferred value. The “true” β and ρ are applied in Eq. (14) to determine the $K_{crit}^2$ values for which T = τ_c. For strong static scattering regions, $K_{crit}^2$ is taken as the minimum of $K_{crit}^2$ and $\beta /\sqrt 2 \; $ to prevent falsely high flow rate measurements due to a near-zero slope produced by very low dynamic scattering (ρ < 0.2; empirically, this limit was only reached with completely static scatters, such as a piece of paper). We run a simple linear interpolation using two points that bound $K_{crit}^2$ to get our final estimate of the correlation time. The first bounding point takes the index, idx, of the last speckle contrast value greater than $K_{crit}^2$ for each pixel. The next index, idx + 1, then acts as the second bounding coordinate for the interpolation:

5. Results

5.1 Performance on simulated data

We applied mm-MESI least-squares fitting (e.g. to Eq. (13), sm-MESI least-squares fitting (e.g. to Eq. (10), REMI, and sREMI to simulated data produced using Eq. (13) with either the traditional exposure times or log-spaced simulated exposure times from Table S1. Parallel processing with 24 workers was utilized only for the least-squares fitting methods. Testing datasets used Eq. (13) to produce 3D matrices where 1/τ_c was varied from 10°−10⁵ s⁻¹ along one axis, ρ was varied from 0.5-1 along the second axis, and the exposure time, T, comprised the third axis (Fig. S3A). β was held as 0.10 while D_MU was set as 0 or 1 to create separate matrices. To test sREMI, a spatial map for 1/τ_c was produced while holding β at 0.10 and spatially randomizing ρ continuously and D_MU as a binary (Fig. S3B). For testing the impact of noise, artificial uniform noise of ±5% of the true value of the resulting speckle contrast value was applied across the 3D matrix. Initial estimates for least-squares curve fitting used values of β=0.12, ρ=0.9, τ_c= 0.001. We additionally set the initial estimate of D_MU as 0.5 for mm-MESI.

For mm-MESI, fitting of the parameters was highly accurate from the shortest exposure time to the second longest exposure time and maintained this accuracy for both scattering models (Fig. 3(A), S4A). β and ρ began to deviate at correlation times below the shortest exposure time while ρ also deviated above the second longest exposure time. The effects of both were mirrored in errors in the estimation of the correlation time. With sm-MESI, estimations of all the parameters were improved over mm-MESI for multiple-ordered scattering when D_MU = 0 (e.g. as we would find in a vessel) but highly deviated for multiple-unordered scattering when D_MU = 1 (e.g. as we would find in a parenchymal region), as expected (Fig. 3(B), S4A). REMI had decent accuracy for correlation times between 10⁻⁴ and 10⁻² s (Fig. 3(C), S4C) with slightly better estimates for multiple-ordered scattering as compared to multiple-unordered scattering. Estimation of β and ρ were accurate over a narrower range of correlation times as compared to mm-MESI, which, in turn, impacted the range over which the correlation time could be reliably estimated. Using evenly-spaced logarithmic exposure times led to significant improvements in REMI-based estimates of each parameter over the “traditionally used” unequally-spaced exposure times, in particular for the multiple-unordered scattering model (Fig. 3(D), S4D). By applying a spatial Gaussian filter (size = 500² pixels, σ = 250 pixels) to spatially smooth β in sREMI, accuracy for estimations of the correlation time significantly improved and the range of accuracy (within 20% of true value) was extended to the shortest exposure time, the accuracy for the estimation of β was expanded across all exposure times, and accuracy for the estimation of ρ was extended to correlation times below the shortest exposure time (Fig. 3(E), S4E). In summary, REMI, and especially sREMI, produce nearly as accurate an estimate of 1/τ_c, β, and ρ as with the full mixed-model least-squares approach, with slightly higher errors in 1/τ_c and β in regions with high flow speeds (with REMI, but not sREMI) and slightly higher errors in 1/τ_c and ρ in regions with low flow speeds (with both REMI and sREMI).

 

Fig. 3. Performance of processing methods on simulated data. Plots of estimated correlation times, τ_c (left), the normalization coefficient, β (middle), and the fraction of dynamic scattering, ρ (right), each normalized to the true value, vs. true correlation times for (A) mm-MESI, (B) sm-MESI, (C) REMI using typical exposure times, (D) REMI using logarithmically spaced exposure times, and (E) sREMI using logarithmically spaced exposure times. Red traces correspond to multiple-ordered scattering data; blue traces correspond to multiple-unordered scattering data. Shaded regions represent standard deviations. Vertical dashed lines represent the bounds set by the exposure times.

Download Full Size | PDF

The mixed-model least-squares fitting method had a slightly higher accuracy and maintained that accuracy over a broader range of correlation times, as compared to REMI, but took 700-1400 × longer to process with 10⁶ pixels and 15 exposures, depending on the unequal or logarithmic spacing of the exposure times (Fig. 4). When varying the image size from 64²−2048² pixels and number of exposures from 10-25 logarithmically-spaced exposure times, the time for least-squares fitting for both models showed a near-linear dependence on the number of pixels with a small influence by the number of exposures, while REMI was equally dependent on the number of pixels and exposure times. With evenly, logarithmically spaced exposure times, and using the compiled filtering convolution function, the exponent for the dependence on the number of exposures was further reduced to ∼0.5. The resulting correlation coefficients between true and estimated τ_c and the respective processing time coefficients for each method are given in Table 1.

 

Fig. 4. Time required for each processing method on 1000 × 1000-pixel images for commonly used nonuniformly spaced exposure times and for uniform, logarithmically spaced exposure times.

Download Full Size | PDF

Table 1. Correlation coefficients between the true τ_c and the estimated τ_c for each method and the processing time required based on the number of pixels and exposure times as ${A\; N}_{{pix}}^{B}{N}_{{exp}}^{C}$.

View Table

5.2 Comparison between methods in a mouse model of photothrombotic stroke

To test the utility of REMI in vivo, a mouse was imaged on a custom setup under anesthesia with respiratory rate and heart rate monitored by a piezoelectric sensor [33] while inducing a photothrombotic occlusion in an arteriole using Rose-Bengal solution and illumination by a green laser with full-width half-max diameter of ∼500 µm at the brain surface (Fig. 5(A)) [34–36]. The imaging setup consisted of a 785-nm laser diode (LD785-SEV300, ThorLabs) which passes through an acousto-optic modulator (AOM; AOMO 3100-125, Gooch & Housego) driven by a radiofrequency (RF) driver (1110AF-AEFO-1.5, Gooch & Housego) to modulate the diffracted laser power. The beam is then expanded for illumination over the sample. An Arduino Uno microcontroller, which communicated with a custom MATLAB program for image acquisition, set the voltage of the RF driver to adjust the amplitude of the diffracted laser beam and gated the pulses of light and camera acquisition using a specialized circuit board. Laser speckle images were collected through a 4 × 0.28 NA objective (XLFluor, Olympus), 792/64-nm bandpass filter (AVR-optics) and a linear polarizer (ThorLabs) and onto a CMOS camera (aca2040-90umNIR, Basler). Due to camera hardware limitations, for exposures less than 50 µs, the camera was triggered for 50 µs while the laser light was pulsed for the shorter exposure time. The TTL signals to the camera, output voltage to the RF driver, piezoelectric sensor signal, and marked trigger signals indicating green light exposure were recorded on a DAQ board (USB-6001, National Instruments) for synchronization. We applied the different processing methods to image data taken from the mouse during repeated photothrombotic strokes. ROI were selected using LSCI analysis of the 4-ms exposures from the dataset before and after induction of the stroke. Correlation times were determined for mixed-model fitting by using Eq. (13), simple-model fitting by using Eq. (10), and LSCI by using Eq. (8) with the 4-ms exposure images, while setting β as the average value from mixed-model fitting.

 

Fig. 5. Processing of speckle images from a mouse model of photothrombotic stroke. (A) Simplified schematic of multi-exposure speckle imaging setup. A microcontroller controls pulse duration and the amplitude of the transmitted laser power through the acousto-optic modulator (AOM) as well as camera exposure times for speckle imaging. The data acquisition (DAQ) board records signals for post-imaging synchronization. Neutral density (ND) filters set the power of the green laser used to produce the photothrombotic stroke, while optical cleanup filters extinguish residual infrared light to prevent interference with the laser speckle imaging. Maps produced by sREMI before (B) and after (C) photothrombotic stroke (green X), showing correlation time (i), normalization coefficient (ii), fraction of dynamic scattering (iii), and scattering model coefficient (iv) with selected ROIs over arterioles (red), venules (blue), and parenchyma (yellow) indicated on the correlation time maps. Scale bar is 200 µm. Greyscale bars are consistent across B and C. (D) Traces of changes in perfusion from indicated ROIs in arterioles (top), venules (middle), and parenchyma (bottom) as measured using mm-MESI, sm-MESI, LSCI, and sREMI from left to right during photothrombotic stroke (green bars) using a 10-point moving average to minimize fluctuations from heart rate. (E) Scatter plot of inverse correlation times from ROIs placed over arteriole (red), venule (blue), and parenchymal (black) regions between LSCI (left), simple-model fitting (middle), and sREMI (right) and mixed-model fitting. (F) Distribution of relative errors from LSCI (left), simple-model fitting (middle), and sREMI (right) measurements with mixed-model fitting as ground truth.

Download Full Size | PDF

To create a vessel occlusion, we used an initial exposure to green laser light of 60 s at an average power of 1 mW, which induced a partial, transient occlusion, then waited some time before inducing a full occlusion by exposure to an increased average green laser power of 2 mW. Average maps before and after photothrombotic occlusions with overlaid ROIs were produced using 100 frames from sREMI outputs (Fig. 5(B), (C)). Movies of 3-point moving average frames produced using sREMI with a spatial Gaussian filter (size = 101² pixels, σ = 50 pixels) on a 5-point moving average of speckle contrast images are shown in Visualization 1 (partial occlusion) and Visualization 2 (subsequent, full occlusion) with whole field-of-view maps of correlation times (τ_c, left), normalization parameter (β, top-right), dynamic scattering fraction (ρ, middle-right), and the scattering model (D_MU, bottom-right). The fraction of dynamic scattering is seen to decrease slightly during the first, mild photothrombotic stroke (Visualization 1) with the second occlusion producing a significant decrease in the fraction of dynamic scattering (Fig. 5(B), C panel iii) and a notable transition between scattering models due to decreased contribution from multiple, ordered scattering after the surface arteriole was fully clotted (Fig. 5(B), C panel iv) (Visualization 2).

All methods were able to capture flow changes during the transient occlusion with 1 mW and following the full occlusion at 2 mW in the targeted arteriole (Fig. 5(D)), but the estimated magnitude of the flow decrease varied by method. After the 1-mW irradiation, there was a noticeable, transient drop in flow in the targeted arteriole that returned to baseline in a couple minutes (Fig. 5(D), top), indicating a partial occlusion that was released upon termination of the light exposure. Mixed-model fitting and sREMI estimated the flow speed decrease in the arteriole to be 17% and 18%, respectively, while simple-model fitting estimated a 15% reduction and LSCI estimated the drop to be ∼41%. With 2 mW of green light exposure, mixed-model fitting and sREMI measured a reduction of 68% and 66%, respectively, in the targeted arteriole from the initial baseline. Simple-model fitting showed a decrease of 60% while LSCI indicated a reduction of 72% from the initial baseline.

To quantify the degree of agreement between methods, we examined the correlation between the inverse correlation times across models, taking the mixed-model least-squares results as ground truth. The simple-model fitting had the overall highest correlation coefficient to the mixed-model fits (R² = 0.9823; Fig. 5(E), left), but there were large relative errors, with a bimodal distribution: 3.9 ± 2.9% and 20 ± 7.4% for the two peaks, respectively (Fig. 5(F), left). In comparison to the mixed-model, the simple-model fitting tended to overestimate flow magnitude in parenchymal regions and in occluded vasculature, leading to the peak at ∼20% error (Fig. 5(F), middle). LSCI estimations of flow tended to be significantly slower than the mixed-model estimates for all ROI types and had the lowest correlation (R² = 0.7648; Fig. 5(E), middle) and largest errors relative to the mixed-model, again with a bimodal distribution: −52 ± 10% and −76 ± 5.6% (Fig. 5(F), middle). sREMI showed a high degree of correlation with the mixed-model least-squares fit in flow changes for all types of ROIs (R² = 0.9785; Fig. 5(E), right), and had the lowest relative error of 2.6 ± 6%, with a normal distribution for this error (Fig. 5(F), right).

6. Discussion

In this paper, we have devised an algorithm that allows for real-time processing of multi-exposure speckle images without the need for least-squares fitting. We validate the efficacy of the algorithm on simulated data, further improve performance using uniform logarithmic exposure times compared to typically used exposure times, and accurately track changes in a mouse during photothrombotic stroke while using mixed-model least-squares fitting as a ground truth.

While REMI is shown to be much quicker than the least-squares fitting approaches, it comes with the drawback of being more influenced by signal noise. As seen in Fig. 1(D), there is prominent salt-and-pepper noise present due to some pixels having higher variability in speckle contrast at the shorter exposure times, resulting in a noisy estimation of β. By applying a spatial filter to β, due to the slow variations expected across the field of view and as we do with sREMI, the high variability at the shortest exposures is eliminated and thus a significant portion of the salt-and-pepper noise is removed (Fig. 1(D)). Additionally, we note an increase in vessel clarity compared to that of least-squares fitting, as evidenced by less unrealistic spiking across the vessel profile. Averaging multiple subsequent REMI or sREMI output images further reduced fluctuations to appear as in Fig. 5(B), C. While some salt-and-pepper noise persists with sREMI, the noise can be remedied with a 5 × 5, or similar, median filter. We note that in the presence of noise during simulations, the performance of sREMI appears very similar in performance to mm-MESI least-squares fitting (Fig. S4A, E), which is further corroborated by the high correlation and low relative error in the estimates of correlation times for the stroke model (Fig. 5(E)-(F), right). The range of correlation times where we find high agreement to the ground truth using sREMI based on simulations also matches well with reported inverse correlation times found for parenchymal tissue and moderately sized blood vessels in mice [8,12], further emphasizing the utility of this method in research applications requiring real-time evaluation.

Unlike least-squares regression methods, REMI allows for estimation of all pixels in an image in a near-real time, or real-time, manner. The time required for REMI to operate scales near-linearly with the number of pixels and exposure times when using non-logarithmically spaced data, such that a 512 × 512 × 15 matrix takes ∼0.25 s to process. A major limitation in speed comes from the nonuniform moving average filter across exposure times applied to the original data as well as the first-order differential using the uncompiled processing in MATLAB. We improved speed by equally spacing the exposure times on a logarithmic scale and applying a lower-level convolution function to filter the data, which reduces the processing time by half to ∼0.12 s for a 512 × 512 × 15 matrix. The processing time difference between sREMI and REMI is near-negligible (Fig. 4, Table 1), influenced primarily by the size of the Gaussian filter for spatially filtering β. As REMI calculations are done across all pixels in parallel, the algorithm is slowed by memory allocations for large matrices associated with larger images, meaning that least-squares fitting methods would eventually take less time to run, though this only begins to occur at unreasonable image sizes on the order of megapixel to gigapixel side lengths. Of course, the image size for nonlinear least-squares regression methods to match the speed of REMI will also change with available parallel processing capabilities. In this study, we used 24 workers in parallel to match the total number of physical cores of the processing machine. Without parallel processing by multiple simultaneous workers, the difference in processing times would be much larger.

Of the parameters fit, the β term is expected to vary the least between samples and across the image field, and can be approximated using images from a static phantom. The impact on processing time by fixing β, however, was negligible. Similar to previous work [17,18,37], we also found that different static surfaces and surface curvature led to slightly different speckle contrast values and, while the laser wavelength and temperature were stabilized, we noticed small fluctuations in the speckle contrast over time likely due to fluctuations in the laser coherence length (data not shown). Based on these findings, we omitted an initial calibration and allowed for β to vary spatially and temporally for all fitting models.

More commonly used single-exposure LSCI is a fast and simple alternative to the quantitative multi-exposure speckle imaging approach, but it comes with the drawback of not being able to account for changes in static scattering, causing significant errors in flow measurements. During the mild photothrombotic stroke, flow reductions were estimated by LSCI to be much larger than with any other method. The errors are largely due to the inability of LSCI to account for changes in the fraction of dynamically scattered light that comes from moving scatterers and not accounting for a shift toward unordered scattering at the site of a vessel occlusion (Fig. 5(B), C panel iii). We find that a relative reduction in the fraction of dynamic scattering by ∼5%, as we saw during the mild stroke, results in a ∼20% underestimation of flow by LSCI. Not allowing the scattering model to change from more ordered to more unordered scattering, which was observed after the full occlusion at the location of the targeted vessel, results in an additional ∼15% reduction in the LSCI flow estimate. Based on Monte-Carlo simulations, the laser speckle signal integrates scattering signals from up to 700 µm below the tissue surface [31]. As the green laser light did not occlude the surface arteriole and flow speed in the arteriole returned to baseline following termination of the green laser light illumination, it is possible that the initial photothrombotic occlusions were restricted to the capillaries just beneath the surface, leading to fewer moving and more static scatterers but leaving the network largely unaffected. The result would be a decrease in the dynamic scattering coefficient while holding the correlation time near constant or only slightly reduced, leading to errors in LSCI flow speed estimation. The sREMI approach applied to this data tracked these changes and produced estimates that closely matched the mixed model, least squares approach.

Machine learning approaches can simplify and speed up computation of complex problems and are generally great for classification purposes and pattern recognition, and have been applied for the estimation of flow from MESI images [19,20]. While versatile, machine learning approaches often require extensive training data and can risk overfitting to the training data, leading to errors when confronted with information outside the training sets [21]. Due to the “black box” behavior of machine learning approaches, it can also be difficult or impossible to determine the origin of any resulting errors [21]. Finally, training the algorithm on data from one system will likely lead to system-dependent accuracy, requiring retraining of the algorithm after modifying the system [21]. As such, having an algorithm with known inputs and processing steps, and that is system-independent can lead to more reliable outputs and easier identification of errors, even if the time required for processing is slightly longer.

7. Conclusion

This paper demonstrates a method for quasi-analytically approximating the correlation times for MESI images in a rapid and reliable manner, with close agreement with the output from least-squares approximations, without the need for pixel-by-pixel curve fitting, removing a significant bottleneck for real-time viewing. The up to 1500 × reduction in processing time allows for monitoring of blood flow changes in real time in clinical and research settings without extensive equipment or machine learning. With proper parallelization, an image set can be recorded, while converting a previous set to speckle contrast images and processing to MESI images then have all outputs displayed with no noticeable delays. All MATLAB code for the implementation of REMI and sREMI, as well as the example data from Fig. 1, is available on GitHub [38].