1. Introduction

A laser speckle is a random interference pattern generated by the interaction between coherent illumination and scattering media [1–3]. Analysis of speckle patterns provides information on the scattering media, such as a spatial distribution of dynamics of the moving scatters like blood cells inside living biological tissue. Various non-invasive speckle imaging methods using array-type photodetectors, such as charge-coupled device (CCD) [4] or complementary metal oxide semiconductor (CMOS) image sensors [5], have been developed. For example, laser speckle contrast imaging (LSCI) is a blood flow imaging technique that measures relative perfusion from temporal fluctuations of speckle patterns [6]. Fluctuations of speckle patterns are converted into speckle contrast, the ratio of standard deviation to mean of speckle intensities in a local region of interest (ROI). Since speckle contrast varies with the speed of motion relative to the integration time of the detector, blood flow distributions can be easily visualized within the imaging field of view [7,8]. However, conventional LSCI lacks absolute flow measurement and underestimates flow rates in the presence of static scattering [9]. Flow rate calibration under multiple scattering events is also unavailable in LSCI.

Several methods based on the temporal analysis of time-varying speckles were developed to resolve speckle dynamics in more complex scattering regimes. For example, temporal LSCI reduces static scattering by calculating speckle contrast from sequential speckle images [10], along with the improvements of complex speckle evolution models [11]. Multi-exposure speckle imaging (MESI) uses speckle contrasts at different camera integration times to quantify the decorrelation time ($\tau _c$) of time-varying speckles under the presence of static scattering [12,13]. MESI has been developed for fast and accurate in vivo blood flow imaging with compact optical implementations [14], quasi-analytic solutions [15], and machine learning approaches [16]. To directly resolve the temporal decorrelation behavior of the speckle fields, dynamic light scattering imaging (DLSI) using high-speed cameras was developed [17]. DLSI directly measures the temporal intensity correlation of speckles at each position in the sequence of speckle images, which can be used to fit more complex scattering models with different dynamics regimes [18]. DLSI brings insight into light scattering and particle dynamics inside biological tissues [19–21].

Recently, we proposed a simple and cost-effective temporal speckle correlation imaging technique named rolling shutter speckle imaging (RSSI) that measures the direct temporal speckle decorrelation behavior in a single shot by taking advantage of the rolling shutter scheme in consumer-grade CMOS image sensors [22]. RSSI can obtain the decorrelation time and the amount of static scattering from a snapshot speckle image by fitting the row-by-row intensity correlations (RICs) of elongated speckle patterns with the spatiotemporal speckle intensity correlation (STIC) model. Thanks to the short row times of modern CMOS image sensors, RSSI can easily demonstrate a sufficient temporal resolution equivalent to the sampling rate of several thousands of frames per second and the extensive range of measurable blood flow rates provided by the microsecond-level row time of the rolling shutter image sensor. RSSI also enables continuous monitoring of sample dynamics at video rates and provides wide-field 2D in vivo blood flow and tissue distributions with high temporal and spatial resolutions using a small amount of image data. However, the original STIC model of RSSI only applies for correlation times longer than the camera exposure duration and demonstrates limited fitting accuracy due to the loss of information. In addition, while the original theoretical model considered the presence of static scattering, the fidelity of the model was not validated, and the accuracy of the estimated scattering parameters has not yet been experimentally verified.

In this paper, we have augmented RSSI experimentally and theoretically and performed a study about its theoretical limit. We first introduce a general STIC model that covers the entire correlation time range, enhancing the fidelity of the model fitting and scattering parameter estimation. We also present simulation results that inform the measurability of RSSI at every imaging condition. Then, we report on the phantom experimental results to verify that the general model can accurately measure the speckle decorrelation times over a wide range of dynamics in the presence of static scattering. With these improvements, we believe that RSSI can serve as a robust, cost-effective imaging technique that resolves static scattering problems and enables quantitative measurement of deep blood flow for biomedical applications.

2. RSSI theory

2.1 Working principles of RSSI

The working principle of RSSI is illustrated in Fig. 1. A rolling shutter CMOS image sensor sequentially records the intensities at every row with a duration of the camera exposure time $T$. The short time interval between adjacent readouts is called row time $t_r$ and has a constant value of 12 $\mathrm {\mu s}$ in our experimental setup. To measure temporal decorrelations of a speckle field on a single image taken by a rolling shutter CMOS image sensor, we need to elongate speckle patterns toward the readout direction (Fig. 1(b)) such that the evolution of each speckle over time can be measured over different rows of the image sensor. In RSSI, we utilize a narrow elliptical aperture (Fig. 1(a)), and its objective numerical aperture (NA) and laser wavelength determine the speckle length. Since speckle fields evolve rapidly within dynamic scattering media as the image sensor records intensities line-by-line, the faster velocity of the moving particles results in the shorter length of the captured speckle patterns in RSSI (Fig. 1(c)). Also, the amount of static scattering inside the media affects the length of the speckle patterns (Fig. 1(c)). Combinations of these two factors determine decorrelation rates of RICs calculated from snapshot speckle images (Fig. 1(d)). To obtain the local velocity of moving scatterers, we should accurately estimate the decorrelation time from regional RICs by extracting the contribution of static scattering. For this, we need a general spatiotemporal intensity correlation model that applies to all correlation time ($\tau$) ranges and scattering levels.

 

Fig. 1. Overview of RSSI. (a) Optical implementation of RSSI composed of a 4-f imaging system including a narrow elliptical aperture on the Fourier plane, a laser diode, tissue and blood-mimicking phantoms, and an elongated point spread function generated on the image plane. (b) Formation of a snapshot speckle image by a rolling shutter image sensor. (c) Snapshot speckle images taken under different velocity and static scattering combinations. (d) Calculated RICs corresponding to the different velocity and static scattering combinations. Fluctuating intensity ratio $\rho$ and decorrelation time $\tau _c$ are obtainable by fitting the RICs to the general STIC model.

Download Full Size | PDF

On a snapshot speckle image, the RIC between the $k^{\mathrm {th}}$ and $\left (k+n\right )^{\mathrm {th}}$ rows can be calculated as:

Here, $\langle \cdot \rangle$ indicates ensemble averaging over the calculation window. The two different rows have the temporal lag $\tau$ and spatial displacement $\upsilon$ as functions of row interval $n$, related by the row time $t_r$ and the pixel size $p$ of the image sensor, respectively (Fig. 2(a)).

 

Fig. 2. Schematic representation of the general spatiotemporal intensity correlation model. (a) The row-by-row intensity readout scheme of a rolling shutter CMOS image sensor. $T$ is the camera exposure time, and $t_r$ is the row time. $\tau$ and $\upsilon$ are the correlation time and spatial distance between two different rows, and both can be represented in terms of the row interval $n$. (b) Schematic and mathematical expressions of the RIC of $I(0)$ and $I(\tau )$. Element-wise intensity correlation in cases of (c) $\tau \ge T$ and (d) $\tau < T$. The connecting lines in the same color demonstrate identical element-wise intensity correlations.

Download Full Size | PDF

The spatiotemporal intensity correlation between the row-by-row speckle measurements has been modeled to obtain the parameters of the speckle dynamics from the RICs [22]. The reported STIC model considers both the temporal correlation of the dynamic speckle field and the spatial correlation given by the shape of the speckles and calculates the intensity correlation between two RSSI measurements that are separated in both time and space. In addition, we have derived the STIC model in consideration of the finite exposure time of the camera for accurate fitting of the experimental data. However, the original model works only for correlation times ($\tau$) longer than the camera exposure time $T$, and using longer $T$ reduces usable RIC data points for estimations of $\tau _c$. In fact, this problem is unique to RSSI, where the integration time can be longer than the sampling interval due to the independent measurements of each row. In RSSI, we can utilize this to improve the signal-to-noise ratio (SNR) of the speckle intensity measurements or increase the measurement range of the speed of the dynamic scattering media. However, due to the lack of the STIC model for when $\tau <T$, the initial data points of the RIC, where the spatial correlation is still high and the effect of the temporal decorrelation is dominant, had to be discarded. This restriction results in low fitting accuracy, especially for rapidly evolving speckle fields with small $\tau _c$ and mixed speckle fields with low dynamic fraction, and also reduces robustness in measurements requiring long $T$, such as slower flow rates and low SNR measurements such as deep blood flow imaging. To resolve this issue, we derive the STIC model for all correlation time ranges to make RSSI applicable for blood flow imaging in various environments.

2.2 General spatiotemporal intensity correlation (STIC) model

To resolve the limitations of the previous STIC model and improve the fitting accuracy, we derive the intensity correlation functions for all correlation time ranges, covering both $\tau \ge T$ and $\tau <T$. In the two measurements of $I(k)$ and $I(k+n)$, the speckle field decays due to the temporal displacement by $\tau$ and spatial displacement by $\upsilon$. The RIC corresponds to the intensity correlation as a function of spatial and temporal displacement, $g_2(\tau, \upsilon )$. To consider the effect of finite exposure times, we consider each intensity as a summation of $N$ consecutive instantaneous intensity components (Fig. 2(b)). If $N$ is large enough, the intensity components and the time gap between adjacent intensity components are infinitesimal, and intensities $I(k)$ and $I(k+n)$ can be expressed as a temporal integration of $I$ within the time ranges of 0 to $T$ and $\tau$ to $\tau +T$, respectively. The multiplication of $I(k)$ and $I(k+n)$ would correspond to the sums of the multiplications between each instantaneous intensity element, as shown in Fig. 2(b). The spatial intensity correlation only depends on the row difference, which is constant for all pairs of instantaneous intensity elements. Assuming the speckle ergodicity [1], the temporal correlation only depends on the absolute time difference between the instantaneous intensity pair. When $\tau \ge T$, the time lags between two instantaneous intensity elements can be written as $\tau +(j-i) \frac {T}{N}$, where $i$ and $j$ are the indices of instantaneous intensity elements (Fig. 2(c)) within $k^{\mathrm {th}}$ and $\left (k+n\right )^{\mathrm {th}}$ rows, respectively. We get multiple sets of instantaneous $g_2$’s with the same time lags ranging between $\tau -T$ to $\tau +T$. By extending these summations into a continuous integral, the numerator of the temporal intensity correlation function for $\tau \ge T$ can be written in terms of the instantaneous $g_2(\tau, \upsilon )$ as follows:

The instantaneous $g_2(\tau )$ is related to the field correlation $g_1(\tau )$ via the Siegert relation [2,3,23]. Here, we assume cross-spectral purity [24], and the field correlation of speckles can be expressed as a product of the temporal and spatial correlation. The speckle field dynamics is typically modeled as a negative exponential decay with the characteristic speckle decorrelation time of $\tau _c$. The spatial correlation between the two measurements decays in the form of the normalized Bessel function of the first kind, which is the transfer function of the elliptical aperture. As a result, the spatiotemporal field correlation can be denoted as:

With the Siegert relation, the instantaneous $g_2$ can be expressed in terms of spatiotemporal correlation terms. Integrating the Eq. (2), the intensity correlation function at the detector integration time $T$ for the correlation times $\tau \ge T$ can be written as:

On the other hand, when $\tau <T$, the time difference between the instantaneous intensities, $\tau +(j-i) \frac {T}{N}$, can be positive or negative. When the time difference is negative (e.g. blue lines in Fig. 2(d)), the intensity correlation can be written as $g_2 \left ( (j-i)\frac {T}{N} -\tau \right )$. By considering the time ranges between $\tau$ to $T$ where the time differences are negative, the finite summations can be expanded to continuous integration of the intensity correlations as follows:

Integrating the equation above, the normalized temporal intensity correlation function of dynamic speckle fields for $\tau <T$ can be written as:

When substituting $\tau =T$ to Eqs. (4) and (6), we can observe that the results are identical and the intensity correlation function is continuous. By using the relationships of correlation time $\tau$ and spatial displacement $\upsilon$ with row interval $n$ (Fig. 2(a)), the spatiotemporal intensity correlation $g_2(\tau,\upsilon )$ can be written as a row correlation function of $g_r(n)$.

2.3 STIC model in the presence of static scattering

Then, to derive the $g_r(n)$ model with finite exposure time with the inclusion of static scattering, we first need to consider the correlation model of the mixed speckle fields [26]. The temporal correlation of the mixed speckle field can be expressed as a weighted sum of dynamic-only and static speckles:

When $\rho = 1$, the above equations yield the same results as Eqs. (2) and (6). With the new STIC model, we can utilize all RIC data points for the curve fitting, which can more reliably estimate the speckle decorrelation time and the dynamic scattering fraction. The experimental results using tissue and blood-mimicking phantoms described in Section 4.2 verify the validity of the general model and the robustness of RSSI in measuring blood flows in the presence of static scattering.

2.4 Relationship between RSSI and speckle contrast $K$

The STIC of the dynamic speckle field measured in RSSI includes the information about the speckle contrast $K$. It is known that the speckle contrast is related to the intensity autocorrelation of the speckle [9]. From Eq. (8), the intensity autocorrelation at $\tau = 0$ yields the following:

3. Measurement range of RSSI

Based on the newly derived STIC model, we explore the parameter space of RSSI to understand the measurable range of speckle dynamics. The range of measurement of RSSI depends on its ability to estimate the correct values of $\tau _c$ and $\rho$ from the RICs. The decorrelation behavior of each RIC is determined by multiple parameters, including $\tau _c$, $\rho$, $T$, and $\beta _s$. Speckle fields with different sample dynamics show distinct RIC curves; however, the sensitivity of the RIC curve to the speckle parameters may vary depending on the parameter values, especially $\tau _c$. For example, when $\tau _c$ is relatively small, i.e., in the order of the row time of the sensor, RIC drops rapidly in the first few rows, and the decorrelation curve is highly sensitive to the small changes in $\tau _c$. On the contrary, when $\tau _c$ is large, exceeding the time corresponding to the speckle length, the RIC is less sensitive to changes in $\tau _c$, and curve-fitted values $\tau _c$ are more prone to errors. The range of $\tau _c$ that can be reliably estimated varies significantly with the camera integration time and $\rho$, as well as with the detector row time and the speckle length. To quantify these ranges, we explore the parameter space of the RSSI by computing the theoretical model with various combinations of parameters to assess the measurement range and errors of $\tau _c$ and $\rho$.

We study the behavior of RIC curves by computing and plotting Eq. (8) for various combinations of $\rho$, $\tau _c$ and $T$. We have fixed the speckle length to 68.45 $\mathrm {\mu m}$, which determines the width of the sombrero function determined by the vertical NA of 0.0115 used in all of our experiments. We also fixed the $\beta _s$ to 0.7, an experimental value we obtained from the speckle imaging experiments using the laser diode in our current setup (LD785-SEV300, Thorlabs, Inc.). Then, the maximum value of $g_r(0) = 1 + \beta _s = 1.7$, and the corresponding maximum $K$ is 0.83 according to Eq. (10). For a specific $g_r(0)$, there exist numerous ($\rho$, $\tau _c$) pairs that satisfy Eq. (9) when $\beta _s$ and $T$ values are determined. We first set the range of $\rho$ from 0.01 to 0.99 and then calculated $\tau _c$ values corresponding to each $\rho$ by solving Eq. (9) using the Vpasolve function in MATLAB. By substituting the obtained ($\rho$, $\tau _c$) pairs to Eq. (8), we have calculated and plotted many different RICs that have the same $g_r(0)$ value. Each set of curves indicates different speckle dynamics, which would not have been distinguished by single spatial contrast, $K$. However, RSSI can differentiate from differences in the temporal decorrelation behavior.

Figure 3(a) shows sets of spatiotemporal intensity correlation curves grouped with the same $g_r(0)$ (or speckle contrast, $K$), where $\rho$ varies from 0.1 to 0.9. The exposure time was set to $T = 12$ $\mathrm {\mu s}$, which equals the row time and the minimum exposure time of our image sensor. Within the curves in the same set, a larger $\rho$ corresponds to a higher dynamic scattering fraction and thus decays faster to 1. The values of $\tau _c$ corresponding to each combination of $\rho$ and $g_r(0)$ are plotted in Fig. 3(b). We observe that for a fixed $g_r(0)$, or a fixed speckle contrast, $\rho$ and $\tau _c$ are approximately linearly related. In the upper region of Fig. 3(b), which corresponds to a high speckle contrast $K$, $\tau _c$ shows a steep increase along with increasing $\rho$, covering a wide range of $\tau _c$, indicating relatively slower sample dynamics. However, the corresponding STIC curves in the rightmost part of Fig. 3(a) reveal that the differences between the STIC curves are very small, indicating that curve fitting may be difficult and result in a large measurement error of $\tau _c$. On the other hand, the curves on the lower left part of Fig. 3(a) show a wide spread of STIC curves, which can be easily distinguishable. However, the corresponding $\tau _c$ range in the lower region of Fig. 3(b) is very narrow, and the differences are minimal.

 

Fig. 3. The measurement range of RSSI. (a) STIC curves with $T = 12$ $\mathrm {\mu s}$ calculated from Eq. (8). The curves are grouped with same $g_r(0)$ values from 1.60 to 1.69 with $\rho$ from 0.1 to 0.9 (b) $\tau _c$ values corresponding to every $\rho$ and $g_r(0)$ values at $T = 12$ $\mathrm {\mu s}$. (c) STIC curves with $T = 36$ $\mathrm {\mu s}$ with the same $g_r(0)$ and $\rho$ as (a). (d) $\tau _c$ values corresponding to the $\rho$ and $g_r(0)$ values for $T = 36$ $\mathrm {\mu s}$. Colorbars in (b) and (d) are shown in the log scale.

Download Full Size | PDF

Applying a longer exposure time of $T = 36$ $\mathrm {\mu s}$ will increase the upper limit of $\tau _c$ range to be measured, as speckle intensity correlation in Eq. (8) is a function of $\frac {\tau _c}{T}$. As expected, the corresponding $\tau _c$ values are approximately three times larger in Fig. 3(d) than in Fig. 3(b). However, increased exposure time reduces RSSI’s ability to resolve the temporal decorrelation behavior of the speckle field. Accordingly, the gap between the STIC curves is reduced for the same $K$ in Fig. 3(c) compared to Fig. 3(a). Further increasing $T$ will result in a small increase in the upper limit of $\tau _c$ at the expense of the RIC resolvability and fitting accuracy. It should be noted that these upper limits in $\tau _c$ are highly dependent on the static speckle length, determined by the size of the elliptical aperture used in the imaging system. Note that all curves and values plotted in Fig. 3 are based on the static speckle length of 37 pixels, corresponding to 444 $\mathrm {\mu s}$ in our detector. If we increase the length of the speckle, the upper limit of the measurable $\tau _c$ will also increase.

4. Experimental results

4.1 Experimental setup

We have implemented the RSSI setup with tissue and blood-mimicking phantom preparations to quantify ($\rho$, $\tau _c$) values under various conditions. For imaging, we have used a USB camera (acA4024-29um, Basler AG) with a 12-megapixel monochrome rolling shutter CMOS image sensor that has a pixel size of 1.85 $\mathrm {\mu }$m and a constant row time of 12 $\mathrm {\mu s}$. The 4-f imaging system contains a $f = 50$ mm tube lens (C VIS-NIR Series, Edmund Optics Ltd.) mounted to the camera and an $f = 50$ mm objective lens (AC254-050-B-ML, Thorlabs, Inc.). We fabricated an elliptical aperture corresponding to a speckle length of 30 pixels in the sensor readout direction and placed it on the Fourier plane of the 4-f imaging system (Fig. 1(a)). To improve the image contrast, we mounted a film polarizer (LPNIRE100-B, Thorlabs, Inc.) between the elliptical aperture and the objective lens. The coherent 785 nm light from a volume holographic grating-stabilized laser diode (LD785-SEV300, Thorlabs, Inc.) was illuminated to the sample plane with a 200 $\mathrm {mW/cm}^2$ power density.

To implement blood flow under human skin, we fabricated tissue-mimicking phantoms with polydimethylsiloxane (PDMS) and titanium dioxide ($\mathrm {TiO_2}$) nanopowders (Sigma-Aldrich, Inc.), according to the recipe published in Ref. [27], which produces a reduced scattering coefficient ($\mathrm {\mu _s}$’) of 2 $\mathrm {mm}^{-1}$. We generated a mixture by adding 85 g of PDMS base, 1 g of $\mathrm {TiO_2}$ mixture, and 20 mL of Toluene [28]. The mixing was done for three days using magnetic stirrers until the Toluene was fully evaporated. Then, to make the size of the $\mathrm {TiO_2}$ nanoparticles homogeneous, we placed the mixture in a sonicator for another three days. Finally, we added 15 g of curing agent to fabricate the complete 1% PDMS/$\mathrm {TiO_2}$ solution, which was cured in tissue phantoms of different thicknesses between 0.12 and 1.56 mm. All phantoms were baked in a heated oven at 75 $^{\circ }$C for 12 hours.

To mimic blood flow with constant flow velocity, we continuously pumped Intralipid 20% solution (I141, Sigma-Aldrich, Inc.) through a single-channel flow cuvette with a cross-sectional area of 20 $\mathrm {mm}^2$ (FCTH13, Cuvet, Co.) using a syringe pump (NE-4000, New Era Instruments). By choosing the PDMS/$\mathrm {TiO_2}$ phantom thickness and controlling the flow rate of the Intralipid solution, we implemented numerous tissue phantom thicknesses and blood flow rates corresponding to distinct $\rho$ and $\tau _c$ values, respectively.

4.2 Tissue and flow phantom measurement results

The main goal of the phantom experiment is to validate the STIC model derived in Section 2.3 with experimental data and to verify the RSSI’s ability to quantify blood flow rates in the presence of static scattering. For robust measurement of sample dynamics, the STIC model should be able to decouple the static and dynamic scattering and correctly estimate both $\tau _c$ and $\rho$ in the experimental data within the measurement range. In the phantom experiment, we controlled the level of static scattering by placing the tissue PDMS/$\mathrm {TiO_2}$ phantoms with 14 different thicknesses from 0.12 mm to 1.56 mm on the flow channel (Fig. 4(a)). At the same time, we provided Intralipid flow rates from 100 to 1600 mL/hr with an increment of 100 mL/hr, corresponding to the linear velocity of 0.14 cm/s to 2.22 cm/s.

 

Fig. 4. Snapshot speckle images and RICs generated from different tissue phantom thickness and flow rate conditions. (a) A tissue phantom on the blood-mimicking phantom flowing through the channel. (b) Snapshot speckle images of tissue and blood-mimicking phantoms. Tissue phantom thicknesses are 0.12, 0.50, and 0.92 mm at each row from top to bottom. Flow rates are 1600, 800, and 200 mL/hr from left to right at each column. (c) RICs (closed circles) and fitted STIC curves (dotted lines) of the snapshot speckle images (1) to (9). With $t_r = 12$ $\mathrm {\mu s}$ and $T = 36$ $\mathrm {\mu s}$, first three data points $\left ( n=0,1,2 \right )$ correspond to the range $\tau <T$ and the rest of the data points correspond to the range $\tau \ge T$ of Eq. (8) in each plot. All scale bars correspond to 100 $\mathrm {\mu m}$.

Download Full Size | PDF

For each tissue phantom thickness and flow rate combination, we captured RSSI snapshot speckle images over $1,000 \times 1,000$ pixels with $T = 36$ $\mathrm {\mu s}$. The imaging system’s focus is maintained on the surface of the tissue phantom layer, and we captured sequences of images at each condition to account for unwanted measurement errors or fluctuations in the exact flow velocity. Figure 4(b) shows snapshot speckle images with different tissue phantom thicknesses (0.12, 0.50, and 0.92 mm in rows) and flow rates (1600, 800, and 200 mL/hr in columns). When comparing the speckle images with the same tissue phantom thickness, we can observe that the speckle length decreases with a faster flow rate. Also, with the same flow velocity, a thicker static scattering layer results in longer overall speckle lengths. On the other hand, the horizontal lengths of the speckles remain relatively constant over all the measurements, with full width at half maximum (FWHM) between 2.20 to 2.35 pixels, with small variations due to a fine drop in the speckle contrast.

We used the same process as Ref. [22] to analyze the speckle images. On a single speckle image, we calculated the regional RIC within the ROI of $50 \times 50$ pixels and performed the ROI scan throughout the entire image with a step size of 10 pixels. Since the sample is spatially uniform, we obtain the representative RIC by averaging the RIC curves within the imaging field of view. Figure 4(c) shows the RICs obtained from the snapshot speckle images in Fig. 4(b) (closed circle). The RICs clearly visualize the differences in the speckle dynamics between samples with different flow rates and/or dynamic scattering fractions. Then, we fitted each RIC with the STIC model, as shown in the dotted line in Fig. 4(c). The improved STIC model in Eq. (8) can utilize all data points in the RIC and demonstrates an excellent fit accuracy of $R^2 = 0.999$ in all measurements, regardless of the thickness of the phantom tissue and the flow rate. We believe that this fitting accuracy shows the fidelity of our theoretical model, especially for data points with $\tau <T$. In addition, we can observe increasing values of $g_r(0)$ for samples with slower flow rates or larger tissue phantom thicknesses, which is consistent with the expected behavior of the speckle contrast $K$ as a function of sample dynamics and the amount of static scattering [26].

Figure 5 shows our phantom experiments’ $\tau _c$ and $\rho$ estimation results. We used 50 raw speckle measurements under the same condition, from which the average and standard deviation of the $\tau _c$ and $\rho$ values were obtained. Figures 5(a) and 5(b) plot the fitted $\tau _c$ and $\rho$’s for selected tissue phantom thicknesses and flow rates, respectively. As a reference, we also measured the speckle dynamics of the flow phantom without the static tissue phantom, which was fitted with the dynamic-only speckle model with $\rho = 1$. In Fig. 5(a), the estimated $\tau _c$ for all tissue phantom thicknesses is inversely proportional to the flow rate, as expected. The trend is consistently observed regardless of the tissue phantom thickness for the flow rates above 300 mL/hr. For data points below the 200 mL/hr flow rate, the decorrelation times were underestimated for most measurements, except for the case with large phantom thicknesses where $\tau _c$ were highly overestimated. This failure at slow flow rates is expected due to the limited RIC resolvability within the high $K$ region, as shown in Fig. 3(c). Also, the lower limit of the speckle decorrelation time that can be measured with RSSI is a function of $t_r$ and the speckle length. In this experiment, our speckle length, as defined by Abbe’s criterion, was 37 pixels, which corresponds to 444 $\mathrm {\mu s}$. In fact, the FWHM of the autocorrelation of the static speckles is roughly 18 pixels, or 216 $\mathrm {\mu s}$; thus, any decorrelation that happens at a longer time scale than these measurement lengths cannot be resolved with RSSI. On the other hand, at sufficiently high flow rates, small differences in the flow velocity are clearly resolved with all phantom thicknesses up to 1.56 mm. However, we note that the values of $\tau _c$ were relatively underestimated for thick phantoms, indicating faster sample dynamics. We believe this may be due to the differences in the dynamic scattering regimes across the samples.

 

Fig. 5. Decorrelation time $\tau _c$ and fluctuating intensity ratio $\rho$ estimated from tissue and blood-mimicking phantoms with different thicknesses and flow rates. (a) The estimated $\tau _c$ values vs. flow rate plots of several tissue phantom thicknesses. (b) The estimated $\rho$ values vs. tissue phantom thickness of several flow rates. Every $\rho$ and $\tau _c$ value has been averaged from 50 consecutive snapshot images. 3D plots of the $\rho$ and $\tau _c$ values estimated from (c) the previous STIC model from Ref. [22] and (d) the new STIC model, corresponding to all 14 tissue phantom thicknesses and 16 flow rates.

Download Full Size | PDF

The fitting results for $\rho$ are shown in Fig. 5(b). $\rho$ values consistently decrease as the tissue phantom thickens for the flow rates $\ge$ 300 mL/hr. On the other hand, for the slow flow rate of 100 mL/hr, $\rho$ is incorrectly estimated for all thicknesses, consistent with the failure of $\tau _c$ estimation at this flow rate in Fig. 5(a), showing the lower limit of flow rate measurement for the current RSSI settings. Also, at 200 mL/hr, $\rho$ is overestimated for thick tissue phantoms. We also note that above the thickness of the tissue of 1.5 mm, the estimated values of $\rho$ start to fail at higher flow rates, which may indicate that these measurements are close to the limit of our detection depth or the lower limit of detectable $\rho$. With our current measurement window size of 50 pixels, we compute the statistics of approximately 25 speckles, indicating that the theoretical limit of the measurable depth will be the thickness where at least one of the speckles within the window is dynamically scattered. Practically, the accuracy of the fit will be low when $\rho$ is small, reducing the maximum detection depth.

Figures 5(c) and 5(d) show the fitting results of our experimental data with the previous STIC model [22] and the new STIC model in Eq. (8), respectively. With the previous model, the loss of the initial data points in the RIC results in a reduced range of reliable fitting. In Fig. 5(c), RICs with large phantom thicknesses with slow flow rates were not properly fitted due to the reduced resolvability of the speckle dynamics. However, in Fig. 5(d), the measurement range of $\rho$ and $\tau _c$ was extended, and most of the data points were well matched and showed good agreement with the phantom thicknesses and flow rates. Furthermore, $\rho$ and $\tau _c$ are no longer overestimated when using the new STIC model. We believe that including data points below the exposure time improved the accuracy of the fitting. Also, the newly derived STIC model is more physically accurate than the estimated model in Ref. [22].

Finally, we verify the validity of the new STIC model using the relationship between the speckle contrast $K$ and the RIC. Figure 6 shows the comparison of $K$ for all experimental data, measured from (a) speckle statistics, (b) RIC, and (c) from the STIC model with the estimated $\rho$ and $\tau _c$ values. We first calculated the speckle contrast $K$ for every tissue thickness-flow rate pair by applying the fundamental definition of $K$ to every 50-by-50-pixel window of the speckle image and averaging all ROIs (Fig. 6(a)). We can also calculate $K$ from RICs using the derived relationship between $g_r(0)$ and $K$ in Eq. (10), as shown in Fig. 6(b). We can also calculate the values of $K$ from the estimated values of $\rho$ and $\tau _c$ from model fitting (Fig. 6(c)). Good agreement between the three measurements of $K$ verifies the STIC model in Eq. (8) and Eq. (10). The spatial contrast was calculated from a single image in each experiment, while $g_r(0)$ and $K(\rho, \tau _c)$ are obtained from 50 images, thus showing a smoother profile. We also note that the $K$ estimations agree well even in regions with incorrect $\rho$ and $\tau _c$, for example, with slow flow rates, indicating that $\rho$ and $\tau _c$ are always estimated in pairs accordingly by the initial values of RIC, which correspond to $K$.

 

Fig. 6. Comparison of $K$ values obtained under different tissue phantom thicknesses and flow rates. (a) $K$ values obtained from $K=\frac {\sigma }{\left\langle I\right\rangle}$. (b) $K$ values obtained by substituting $g_r(0)$ values of RICs into Eq. (10). (c) $K$ values calculated by substituting the fitted $\rho$ and $\tau _c$ values into Eq. (9).

Download Full Size | PDF

5. Conclusion and discussion

In this paper, we have derived and verified an improved speckle spatiotemporal intensity correlation model for RSSI, which provides robust $\rho$ and $\tau _c$ estimations for experimental data over a wide range. The STIC model for the first few data points ($\tau <T$) in RSSI allows us to utilize the initial decorrelation behavior of the dynamic speckle field, which improves the robustness of the measurement and prevents overestimation of $\rho$ and $\tau _c$. With the new model, we explored the parameter space of RSSI to assess the measurable range of $\tau _c$ and study the effect of exposure time on the RSSI measurements. We also experimentally validated the STIC model based on phantom experiments with different dynamic scattering fractions and flow velocities. Compared to the previous STIC model, $\rho$ and $\tau _c$ can be measured over a broader range of values, showing the lower limit of $\rho$ near 0.2 and the upper limit of $\tau _c$ near the FWHM of the speckle length. We also validated the model by showing that the speckle contrast can be easily obtained from the RIC.

The simulation and experimental results demonstrate the advantages and limitations of RSSI in dynamic scattering measurements. Thanks to the direct temporal measurement enabled by the rolling shutter, we can discern the static and dynamic scatterings and measure the dynamics when the fluctuating intensity only takes a small fraction of the total intensity, as low as 20%. RSSI’s temporal measurement range is predominantly determined by the row time of the sensor and the length of speckles. With current experimental settings, it is more apt for measuring fast dynamics. Slower dynamics can be measured by utilizing an image sensor with a longer row time or longer speckle lengths at the expense of spatial resolution. Increasing the exposure time of the detector may push the upper limit of $\tau _c$, but perhaps with more error from model fitting because of the reduced spread between the STIC curves.

This paper provides a comprehensive model of RSSI and an understanding of the parameter space of RSSI. The result of this study can serve as a guide to selecting imaging parameters when using RSSI for various applications of dynamic scattering imaging and quantification. Along with the 2D imaging capabilities of RSSI, we believe we can extend our results to improve and optimize RSSI for in vivo blood flow imaging for single-shot measurement with superior performance at measuring subsurface flow. Furthermore, the STIC model under different scattering regimes and their combinations can help with the in vivo measurements. Currently, the model has not yet been derived for the scattering regimes other than $z=1$. We consider deriving the analytic form of the STIC models applicable for RSSI under ordered flow ($z = 2$) and multiple scattering ($z = 0.5$) regimes as potential future work. In addition, the relationship between $\rho$ and the volumetric fraction of blood within the scattering volume also needs to be studied for physiological translation of quantitative blood flow measurement such as speckle plethysmography [29]. These accurate forward models for RSSI can help utilize physics-based learning to estimate better the scattering parameters for the real-time processing of RSSI images. With these future directions, we expect RSSI to be a cost-effective tool for real-time perfusion imaging and a monitoring tool for cardiovascular diseases.