1. Introduction

Blood flow and tissue perfusion are important physiological and pathological parameters in clinical diagnostics. Laser speckle contrast imaging (LSCI) [1] provides full-field, high-resolution, real-time imaging of blood flow and tissue perfusion. As a label-free imaging modality, LSCI measures the relative blood flow based on the laser speckle contrast analysis (LASCA) of backscattered speckle images [2]. Over the last two decades, LSCI has been used to image the cerebral cortex, retina, skin, mesentery, bone joints, and other structures [3–7]. LSCI is usually used in a reflective imaging geometry with full-field illumination in continuous wave (CW) mode [8,9]. More specifically, it uses the coherent light, e.g., a near-infrared laser, to illuminate tissue surfaces, gathers backscattered speckle images and calculates a contrast image that represents the relative blood flow. The contrast value in each pixel is defined as the ratio of the standard deviation to the mean intensity [10], i.e., $K=$ $\sigma / \mu$. This ratio is theoretically related to the blood flow velocity $v$ through Eq. (1) [2]:

Since LSCI utilizes full-field CW illumination, the Brownian motion of scatterers in tissue causes temporal intensity fluctuations in the recorded speckle images [11]. The ordered blood flow in tissue blurs such dynamics and reduces the contrast values. Compared with pointwise or line-scan illumination, full-field illumination causes more multiple-scattering photons to be detected at each pixel.Single- and multiple-backscattered light components thus coexist in speckle patterns. Singly-backscattered light enters the tissue and is backscattered once. Thus, this component primarily penetrates the superficial layer of tissue, with the majority less than half of the transport mean free path ($< l_t/2$) and full range up to $2l_t$. For brain tissue ($l_t \approx 100 \mu m$) [12], $95\%$ reflectively detected signals come from the top $700 \mu m$ ($\sim 7l_t$) [13]. The effective depth of multiple scattering component is thus about $1l_t \sim 7l_t$. The superposition of both components, e.g. the recorded speckle images, always suppresses some of the details in these components.

Biological tissues are complex forms of random media. Coherent waves propagating through a random medium undergo multiple scattering processes with interference phenomena [14]. For transmissive imaging, optically thick tissue is turbid or opaque because of strong multiple scattering [15]. Even in semitransparent tissues, no clear structures can be observed when the transmitted light is multiply scattered. Thus, for in vivo applications, reflective imaging or detection is preferred. To describe light propagation, a traditional radiative transport equation (RTE)-based model can be used [16], but such a model does not account for the interference phenomena in the speckle images, i.e., the coherent addition of backscattered light with random phasors from different paths. Random matrix theory (RMT) has been successfully used to describe the behavior of CW multiple scattering in random media [17,18], such as light focusing and target detection in turbid and opaque samples using acoustic and light waves [19–21]. We recently established a random matrix (RM)-based description of dynamically-backscattered speckle patterns based on trajectory perturbations due to the Brownian motion of scatterers [22]. A strategy for separating the single- and multiple-scattering components in backscattered light was also developed.

In many cases, separating single- and multiple-scattering components is difficult in reflective imaging. Pulse illumination with time gating (e.g., time of flight) [23], and low-coherence light with coherence gating (e.g., optical coherence tomography (OCT)) [24] utilize the path length difference to filter out multiple scattering components. In addition, a source-detector separation strategy was proposed for filtering out single-scattering components (NIRS) [25]. However, these traditional methods cannot be used with full-field CW illumination. Our RMT-based method facilitates the separation of single- and multiple-scattering components in in vivo tissue imaging applications under full-field CW illumination. In this study, we adopted and optimized the RMT-based LSCI of tissue blood flow and functional changes.

2. Theory

Standard LSCI utilizes full-field illumination with a coherent NIR laser (e.g., $785 \mathrm {~nm}$ in this study) in $\mathrm {CW}$ mode. The backscattered speckle pattern is recorded by a monochromatic camera. The exposure procedure of camera is controlled by setting the exposure time. The gray levels in the recorded image $I\left [n_{1}, n_{2}\right ]$ ($n_{1}=1 \cdots N_{1}, n_{2}=1 \cdots N_{2}$, $N_{1} \times N_{2}$ is the image size) are proportional to the light intensities. There are single- and multiple-scattering contributions to each pixel (entry), i.e., $I\left [n_{1}, n_{2}\right ]=I_{S}\left [n_{1}, n_{2}\right ]+I_{M}\left [n_{1}, n_{2}\right ]$. Since multiple-scattering contributes to the speckle pattern diffusively, the entries in $I_{M}$ follow the Gaussian distribution $\mathbb {N}\left (\mu _{M}, \sigma _{M}^{2}\right )$ with mean value $\mu _{M}$ and variance $\sigma _{M}^{2}$ [26]. As a result, the single-scattering components follow the negative exponential distribution $\mathrm {Exp}(\lambda )$ with mean value $1/ \lambda$ [27].

To obtain the dynamically-backscattered speckle images, we obtain the sequential recordings: $\left \{I_{i} \mid i=1, \ldots, T\right \}$. Standard LSCI uses a monochromatic camera with a low frame rate ($50$ fps) and a short exposure time ($5 \mathrm {~ms}$), the interframe time interval ($\geq 15 \mathrm {~ms}$) is sufficiently long compared with the decorrelation time $\tau _{c}$ ($< 500 \mu s$) [28] and thus ensure the temporal sampling is statistically independent [29]. Brownian motion introduces dynamic fluctuations in the raw recorded speckle images. For the multiple-scattering components, temporally independent sampling follows i.i.d.(independent and identically distributed) realizations from $\mathbb {N}\left (\mu _{M}, \sigma _{M}^{2}\right )$. On the other hand, the single-scattering component produces relatively stable patterns due to the presence of in-phase paths with low-rank characteristics. The blood flow is another source of dynamic effects on both the single- and multiple-scattering components. In the following theoretical derivation, for simplicity, we ignore the blood flow effect at first and revisit its effects in the contrast analysis.

For $T$-frame speckle images $\left \{I_{i} \mid i=1 \cdots T\right \}$, each $I_{i}$ is reshaped to a column vector $h_{i}$, with entries $h_{i}[n] (n=1, \ldots, N, N=N_{1} \times N_{2}$) as the $\mathrm {i}^{\text {th }}$ column in the $N \times T$ temporal intensity random matrix $R$. Because $R=R_{S}+R_{M}$, we can further centralize $R$ by $\hat {R}=R-\bar {R}$, where each entry in $\bar {R}$ is the mean value along the same row of $R$. We thus have $\hat {R}=\hat {R}_{S}+\hat {R}_{M}$, where the subscripts $S$ and $M$ denote the single- and multiple-scattering components, respectively.

We next investigated the spectral density of the centralized Wishart random matrix, $E=\hat {R} \hat {R}^{\prime }$. The eigenvalues of $E$ are denoted by $\{s(i)\}$, where $i=1, \ldots, N$ and $s(1) \geq s(2) \geq \cdots \geq s(N)$. The probability density of the eigenvalues can be calculated $\rho (s) \triangleq \frac {1}{N} \sum _{i=1}^{N} \delta (s-s(i))$.

The entries in the centralized multiple-scattering component $\hat {R}_{M}$ follow the i.i.d. Gaussian distribution $\left (\sim \mathbb {N}\left (0, \sigma _{M}^{2}\right )\right )$. The eigenvalue density of the multiple-scattering components $E_{M}$ is well described by the Marčenko-Pastur law with $T \geq N$ and $N,T \rightarrow \infty$ [30]:

When the eigenvalue density of $E_{M}$ has finite $4^{\text {th }}$ moments, its maximum $s_{M}(1)$ [31] and minimum $s_{M}(N)$ [32] eigenvalues converge:

Equation (3) provides the theoretical basis for estimating the intensity variance of the multiple-scattering components using the extreme eigenvalues of $E_{M}$. In practice, only the eigenvalues of $E$ can be calculated, which could be affected by both $E_{M}$ and $E_{S}$. However, the low-rank characteristics of the single-scattering components limit their contribution to the larger eigenvalues of $E$. The eigenvalues of $E$ thus converge to the minimum eigenvalues of $E_{M}$ when the eigenvalues of $E_{S}$ are nearly zero, which was proven by Loubaton and Vallet [33]. We thus can use the minimal eigenvalue of $E$, i.e., $s(N)$, as an estimate for $\sigma _{M}^{2}$ :

Then, we can analyze the trace relation (Eq. (5)) and the sampling variance $\tilde {\sigma }^{2}=\operatorname {tr}\left (\tilde {R} \tilde {R}^{\prime }\right ) / N T$, $\tilde {\sigma }_{S}^{2}=$ $\operatorname {tr}\left (\tilde {R}_{S} \tilde {R}_{S}^{\prime }\right ) / N T$, and $\tilde {\sigma }_{M}^{2}=\operatorname {tr}\left (\tilde {R}_{M} \tilde {R}_{M}^{\prime }\right ) / N T$.

Since the multiple-scattering component $\tilde {R}_{M}$ follows the Marčenko-Pastur law, we have $\tilde {\sigma }_{M}^{2} \rightarrow \sigma _{M}^{2}$ and $\tilde {R}_{M}[n,t] \sim N\left (0, \sigma _{M}^{2}\right ) .$ Then, we can find that the commutative part $\operatorname {tr}\left (\tilde {R}_{S} \tilde {R}_{M}^{\prime }\right ) / M T \sim \mathbb {N} \left (0, \tilde {\sigma }^{2} / N T\right ) \approx 0$. Finally, we can use Eq. (6) to estimate $\sigma _{S}\left (\sigma _{S}^{2} \approx \tilde {\sigma }_{S}^{2}\right )$.

The corresponding mean intensity $\mu _{S}$ can thus be obtained with the exponential distribution of $\mu _{S}=\sigma _{S}$ and the known value of $\sigma _{S}$. Finally, the mean intensity can be obtained: $\mu _{M}=\mu -\mu _{S}$.

The effect of the ordered flow can be modeled by adding a drift component to the original Brownian motion, which blurs the temporally averaged speckle patterns and decreases both $\sigma _{S}$ and $\sigma _{M}$. In standard LSCI, we calculate the contrast values $K=\sigma / \mu$ to quantify the blurring effect of the ordered flow by assuming the heterogeneity of the mean intensities in the full-field imaging area, which can be directly implemented applied to the multiple-scattering contrast images, i.e., $K_{M}=\sigma _{M} / \mu _{M}$. For the single-scattering components, the ordered flow blurs $\sigma _{S}$ and $\mu _{S}$ simultaneously due to the intrinsic exponential distribution. Therefore, the contrast $K_{S}=\sigma _{S}$ represents the relative flow velocity in the single-scattering contrast image. The relations among $K_{S}$, $K_{M}$ and the flow velocity $v$ were validated using both electric field Monte Carlo simulations and in vitro blood flow phantom experiments in the following sections.

3. Methods and experiments

3.1 Electric-field Monte Carlo (EMC) simulation

Radiation transport equations (RTEs) provide a comprehensive description of light propagation in random media. In practice, analytic solutions of RTEs are difficult. As a numerical tool, Monte Carlo simulation has been developed to investigate light propagation [34]. Both RTEs and traditional Monte Carlo methods ignore interference effects in coherent light transportation. To simulate coherent propagation, we need to track both the changes in the light fields of each scattering event and the light trajectories, for example, by using Mie scattering theory as suggested by Xu and colleagues [35]. The interference effects can be modeled as a coherent addition to the electric fields. In this study, we therefore extended the EMC simulation to full-field coherent illumination and the detection of speckle patterns.

To simulate blood flow in biological tissue, we developed a $400\mu m \times 400\mu m \times 400\mu m$ phantom that contained an aqueous solution of homogeneously and randomly distributed Mie scatterers (diameter $a=0.8 \mu m$, volume fraction 50%). A vessel (radius: $r = 70 \mu m$) beneath the surface layer contained the same scatterers but with additional ordered flow motion. A matching boundary was used to simplify the simulation. The illumination light was simulated with a coherent $800 \mathrm {~nm}$ linearly polarized plane wave, which corresponded to a size parameter $x=2.09$, a transport mean free path $l_{t}=2.066 \mu m$, a refractive index $n=1.59$, and an anisotropic coefficient $g=0.685$. In each simulation, we launched $10^{8}$ photons as the full-field incident light and recorded the position and electric field of each scattering event, as well as backscattered light trajectories. To obtain a steady-state CW illumination output, we added all backscattered path fields coherently, resulting in a single speckle image: $I[n_1,n_2] = |\sum E_x[n_1,n_2]+ \sum E_y[n_1,n_2]|^{2}$.

For each scattering event, the local coordinate $\left (\boldsymbol {e}_{1}, \boldsymbol {e}_{2}, \boldsymbol {q}\right )$ was rotated to $\left (\boldsymbol {e}_{1}^{\prime }, \boldsymbol {e}_{2}^{\prime }, \boldsymbol {q}^{\prime }\right )$ by applying the scattering angles $(\theta, \phi )$ (Eq. (7)) , which were sampled from a random phase of the normal distribution $p(\theta, \phi )=F(\theta, \phi ) / \pi Q_{s c a} x^{2}$, where $Q_{sca}$ is the scattering efficiency, $x=ka$ is the size parameter, $a$ is the radius of the particle and $F(\theta, \phi )$ is the scattered light intensity along the direction $(\theta, \phi )$.

Each scattering event altered the incident electric field $\boldsymbol {E}=E_{1} \boldsymbol {e}_{\mathbf {1}}+E_{2} \boldsymbol {e}_{\mathbf {2}}$ to $\boldsymbol {E}^{\prime }=E_{1}^{\prime } \boldsymbol {e}_{1}^{\prime }+E_{2}^{\prime } \boldsymbol {e}_{2}^{\prime }$ through Eq. (8).

To obtain dynamic speckle images with Brownian motion, new trajectories of the same configuration were generated by moving the scatterers based on the random walk model, i.e., $\Delta r(\tau ) \sim \mathbb {N} (0,\sqrt {6D_B\tau })$, where $D_B$ is the diffusion coefficient of the scatterer. Based on the momentum relaxation time and the configurational relaxation time, the time scale used in the simulation was set to $5\times 10^{-4} s$. The backscattered fields were then recalculated, and the speckle image was refreshed in accordance with the above time scale. To simulate the ordered motion, we added constant drifts $v\tau$ to each $\Delta r(\tau )$ with $v \in [1\sim 6] \mathrm {mm} / \mathrm {s}$ inside the flow regions. For each configuration, 10 speckle images were generated and averaged as one recorded speckle image obtained by the camera (exposure time of $5 ms$).

3.2 In vitro phantom flow experiment

We then used an in vitro phantom experiment to test the accuracy and linearity of the flow velocity estimation using LASCA of the single- and multiple-scattering components, respectively. Intralipid (IL) (Kabivitrum Inc., USA) solution the was used as the Mie scattering random medium. The mean diameter of the lipid droplets was $0.7 \mu \mathrm {m}$, and the scattering coefficient $\mu _{\mathrm {s}}$ at $785 \mathrm {~nm}$ was approximately $2 \mathrm {~mm}^{-1}$ at a 2% concentration. The mean free path (MFP) was $0.5 \mathrm {~mm}$ at a 2% concentration. We used a syringe pump to control the flow of the IL solution inside a polyethylene tube (PE-50, outer diameter: $0.97 \mathrm {~mm}$; inner diameter: $0.58 \mathrm {~mm}$ ) at different velocities $(2 \sim 10 \mathrm {~mm} / \mathrm {s})$.

A diode laser ($785 \mathrm {~nm}$, LP785-SF20, Thorlabs, USA) was used as the monochromic coherent light source. The laser beam was expanded by a diffuser to illuminate the surface of the phantoms (Fig. 3(a)). The reflected light was imaged by a monochrome 12-bit CCD camera (SCA640-70fm, Basler Scout, Germany) with a macro lens (AF-S DX Micro $40 \mathrm {~mm}$ f2.8, Nikon). The recorded speckle images were then used to construct the hybrid Wishart RM and estimate the first and second moments of the single- and multiple-scattering components, respectively. The contrast values of the single- and multiple-scattering components were calculated using different window sizes. The averaged contrast value covering the tube area was substituted into Eq. (1) to obtain $\tau _{c}$ ($\beta = 1$). Finally, the estimated velocity $\tilde {v}$ was calculated by $\tilde {v} = a\tau _{c} +b$ and compared with the true velocity $v$. Here $a$ and $b$ are normalization parameters to ensure $\tilde {v}$ in the range of the true velocity.

3.3 In vivo CBF imaging

3.3.1 Steady-state CBF imaging of rats

All animal experimental procedures were performed in accordance with protocols approved by the Animal Care and Use Committee of Shanghai Jiao Tong University. Five adult Wistar rats $(\sim 250 \mathrm {~g}$, female) were used in the current study and randomly assigned into parameter optimization group (n=1), steady-state CBF imaging group (n=2) and functional CBF imaging group (n=2). Rats were anesthetized with sodium pentobarbital (3ml /kg, IP) and placed in a stereotactic frame (David Kopf Instruments, Tujunga, CA, USA). A homeothermic blanket system was used to maintain the rectal temperature of the rats at $37^{\circ } \mathrm {C}$. After a midline incision was made on the scalp, a high-speed dental drill (Fine Science Tools Inc. North Vancouver, Canada) was used to thin an $5 \mathrm {~mm} \times 5 \mathrm {~mm}$ area centered $3.5 \mathrm {~mm}$ lateral to and $3 \mathrm {~mm}$ posterior to bregma.

The in vivo imaging experiment used a well calibrated small animal imaging system in our lab. A 12-bit cooled monochromatic CCD camera (Sensicam SVGA, Cooke, Michigan, USA) with a $60 \mathrm {~mm}$ f/2.8 macro lens (Nikon Inc., Melville, NY, USA) was used to record the laser speckle images $(1280 \times 1024$ pixels, $10 \mathrm {fps}$ ) under He-Ne laser illumination $(632.8 \mathrm {~nm}, 0.5 \mathrm {~mW}$, JDSU, Milpitas, California). In practice, near-infrared wavelength can also be applied to improve the penetration depth. A total of 30 speckle images were recorded. The single- and multiple-scattering components were separated after constructing the Wishart RM. Then, the contrast images of the single- and multiple-scattering components were calculated. The standard temporal LASCA algorithm was also applied to the recorded speckle images to obtain a hybrid contrast image. The results were compared for different temporal and spatial window sizes.

After LSCI, fluorescence imaging was performed to obtain ground truth images of the cortical vasculature. Under anesthesia, $200 \mathrm {~mL}$ rhodamine-dextran tracer (Invitrogen, Carlsbad, California, USA) was injected into the tail vein. Fifteen minutes later, the rhodamine-dextran level in the cerebral cortex was recorded using the same camera but with a bandpass emission filter ($560 \pm 5 nm$) under the fluorescence excitation of a $532 \mathrm {~nm}$ diode-pumped solid-state (DPSS) laser ($5 \mathrm {~mW}$, Nd: YVO4KTP, Beam of Light Technologies, Clackamas, Oregon, USA). To enhance the SNR of the fluorescence signal, we averaged 20 fluorescence emission images to obtain the final fluorescence image.

3.3.2 Functional CBF response to electrical hind paw stimulation

The surgical preparation process was similar to that of the steady-state rat CBF imaging. After thinning the skull area ($5 \mathrm {~mm} \times 5 \mathrm {~mm}$) above the right somatosensory cortex, two needle electrodes (Safelead F.E3.48, Grass. Telefactor, USA) were inserted subcutaneously into the rat’s right hind paw. Each stimulation trial consisted of a baseline period ($5 s$), electrical stimulation ($10 s$, rectangular constant current pulses ($2.5 mA, 0.3 ms, 5 Hz$) and a recovery period ($300 s$). During direct electrical stimulation, the cranial window was illuminated by a He-Ne laser and imaged by the CCD camera to record the functional CBF response. A home-made software is developed to control and synchronize the camera recording and electrical stimulation through USB 3.0 connection. Starting from the baseline period, the images were continuously recorded at $10 fps$ with exposure time of $5 \mathrm {~ms}$. The contrast images for both the single- and multiple-scattering components were analyzed to obtain the CBF response to electrical stimulation.

4. Results

4.1 EMC simulation and in vitro phantom flow experiment

After all (hybrid) trajectories were recorded in the Monte Carlo simulation, we identified the single scattering trajectories as those with single scattering event. The rest are multiple scattering trajectories. Different kinds of trajectories further generate dynamic speckle images. Figure 1(a) shows the temporal LASCA (tLASCA) of the dynamic speckle images of the single-scattering trajectories, multiple-scattering trajectories and hybrid trajectories, in a region of interest that covers the vessel center. tLASCA of all three kinds of trajectories showed good estimates of the simulated velocities within the range of $1 \sim 6 mm/s$, which proved the validity of LSCI of either single- or multiple-scattering components. In addition, estimated velocities from all (hybrid) trajectories demonstrated the highest linearity ($R^{2} = 0.953$) with the velocity, which validated the model consistency of Eq. (1) (Lorentzian velocity distribution [2]) in tLASCA. However, the linearity of estimated velocities from isolated single-scattering trajectories or multiple-scattering trajectories was degraded ($R^{2} = 0.812$ or $0.886$ respectively), which was mainly due to the less consistency of the model and unbalanced sampling of the trajectories ($45\%$ vs. $55\%$).

 

Fig. 1. EMC simulation results. (a) The calculated velocities using the tLASCA method from the single-scattering trajectory, multiple-scattering trajectory and hybrid trajectory, all showing strong linear relationships with the true velocity. (b) The reconstructed velocity using the separation method based on the estimated single- and multiple-scattering components also demonstrated a linear relation. (c) and (d) show the relative estimation errors of the velocity in the multiple-scattering component for different settings of $Q$ and $S$.

Download Full Size | PDF

Figure 1(b) shows the estimated velocities by using only the single- ($R^{2} = 0.931$) and multiple-scattering ($R^{2} = 0.938$) components, separated from the hybrid trajectory with RMT, which results in greater $R^{2}$ values. The RM based method statistically separated the components by estimating the first- and second-order moments. Such a separation worked on all (hybrid) trajectories and thus maintained consistency with the model of Eq. (1). Nonetheless, new models for better linearity should be developed and verified for the separated single and multiple scattering components.

In vitro experiments were performed to further validate the proposed separation strategy. Using a syringe pump, we are able to accurately control the flow (Fig. 2(a)). Figure 2(b) shows the estimated velocity and corresponding ground-truth velocity based on the single- and multiple-scattering components in the in vitro flow experiments. The proposed method demonstrates high linearity and robustness for both single scattering $(R^{2}=0.803)$ and multiple scattering $(R^{2} = 0.696)$ at a window size $S=3$ (in pixel). Compared with the EMC simulation, the in vitro experiment showed slightly less linearity in the estimated velocity due to the line width and polarization limitations of the diode laser, the noise and distortion in the lens imaging, and the velocity spatial distribution in the laminar flow.

 

Fig. 2. (a) The imaging setup used in the in vitro phantom flow experiment. (b) The estimated velocities from the single- and multiple-scattering components vs. the true velocity in the in vitro phantom experiments. (c) Coefficient of determination $R^{2}$ between the estimated velocity and true velocity using different size windows.

Download Full Size | PDF

4.2 Parameter optimization based on EMC simulations and in vitro flow phantom experiments

Theoretically, a larger Wishart RM is more likely to provide a better estimation of the spectral density, which was confirmed by the EMC simulation. In Fig. 1(d), the estimation error for the multiple-scattering component exponentially decreased to less than $1 \%$ as the window size $S$ (in pixel) linearly increased. It is worth noting that linearly increasing the parameter $Q$($Q=T/N$, $N=S \times S$ ) resulted in a damping attenuation of the estimation error for multiple-scattering components (Fig. 1 (c)). Since the temporal resolution is important for in vivo imaging, the range of $2\leq Q\leq 4$ can be applied to achieve both the sufficient accuracy and temporal resolution.

Figure 2(c) shows the effects of the sampling window size on the linearity of the flow velocity estimation in the in vitro flow experiment ($Q$ is always greater than or equal to 3). An increase in the window size ($S=3 \sim 9$, $N=9 \sim 81$) slightly improved the linearity of the flow velocity estimation for the single-scattering component. The linearity of the flow velocity estimation for the multiple-scattering component showed a significant drop at $S=7$ and decreased to $0.413$ at $S=9$. Thus, a large $N$ is not recommended for realistic applications, especially for the reconstruction of the multiple-scattering components.

Both the EMC simulations and in vitro flow experiments suggest an optimized parameter range of $9 \leq N \leq 49$ and $Q \geq 3$. For in vivo flow imaging applications, due to the spatial diversity of biological tissue, a larger spatial sampling window size ($S$) will cause the spectrum of the Wishart RM to be contaminated with inhomogeneous noise, thereby reducing the spatial resolution of the reconstructed contrast image. Thus, the parameter $S$ should be constrained to a limited area. Additionally, because blood flow constantly changes, it is difficult to keep the scattering signal ergodic for longer $T$(stationary flow assumption). Therefore, once the parameter $N$ has been determined, the parameter $Q$ cannot be arbitrarily large.

4.3 Parameter optimization for in vivo blood flow imaging

In the steady-state in vivo blood flow imaging experiment, we obtained optimized parameters for $N$ and $Q$. Figure 3(a) and (h) demonstrate the reconstructed single- and multiple-scattering contrast images for different window sizes $S$ and different numbers of frames $T$. For single-scattering contrast images (Fig. 3(a)), the noise decreased with increasing $S$ (or $T$) when $T$ (or $S$) was fixed. Similar effects were observed in the multiple-scattering contrast images (Fig. 3(h)) when $S<=5$ and $T<= 60$. However, further increasing $S$ and $T$ (see last image in Fig. 3(h)) introduced more noise. Therefore, we proposed a quantitative measure of the noise level, defined as the sample standard deviation of the contrast values in the region of interest (ROI) covering the larger vessel within the imaging window. Each value was then normalized by the mean contrast value of that ROI.

 

Fig. 3. (a) Enlarged blue circle areas reconstructed from the single-scattering component using different values of $T$ and $S$. (h) The enlarged red circle areas reconstructed from the multiple-scattering component using different values of $T$ and $S$. (b)(e) The contrast values along the cross section of the vessel indicated by the white arrow in (a). (c)(f) The contrast values along the cross section of the vessel indicated by the white arrow in (h). The noise level of the contrast values in the vessel regions calculated by the conventional tLASCA method and reconstructed from the single- and multiple- scattering components by applying different values of $S$ (d) and $T$ (g).

Download Full Size | PDF

The relative noise levels ($mean\pm SD$) in the vessel ROIs of the hybrid, single-scattering and multiple-scattering contrast images are shown Fig. 3. Figure 3(d) showed the relative noise levels using different spatial window size $s$. Multiple-scattering contrast image with $S=3$ or $S=5$ had less noise levels ($<0.1$) than that of tLASCA algorithm ($S=1$) from the hybrid speckle images. For the single-scattering image, a larger $S$ decreased the noise level (Fig. 3(d)). For the multiple-scattering image, a window size that was too large (i.e., $S=7$) significantly increased the noise level due to the increased variance in estimating the second-order moment of the multiple-scattering components. In general, a larger spatial window always reduced the spatial resolution and introduced a distorted velocity distribution, making it difficult to detect smaller vessels. Figure 3 (b) and (c) demonstrate the normalized contrast curves crossing the vessels indicated by the arrows in Fig. 3 (a) and (h). The U-shaped pattern of the contrast values at $S=3$, which represent the inverse distribution of the velocity, disappeared at $S=7$. Thus, the window size $S$ was set to $S=3$ ($N=9$) in the following in vivo experiment for optimal imaging of the vessels.

With the fixed $N=9$, the $Q$ value is determined only by the value of $T$ used to construct the Wishart RM. Figure 3 (g) shows the noise level of the hybrid, single-scattering and multiple-scattering contrast images with different values of $T$ and $S=3$. The single-scattering component had a relatively stable noise level that was always greater than that of the hybrid contrast image. The multiple-scattering component had a high noise level when T is relatively small (e.g., $T=15$ and $Q = 15 / 9$) due to the high variance in the spectral estimation of the Wishart RM as shown the first panel in Fig. 3 (h). An insufficient number of frames $T$ also introduced distortions in the normalized contrast value curves (Fig. 3(e) and (f)) along the cross section of the selected vessels. The noise level in the multiple-scattering contrast images could be lower than that in the hybrid contrast images if $T$ was sufficiently large (e.g., $T \geq 30$). It is thus critical to optimize the parameter $Q$ with minimal $T$ in functional CBF imaging. A large $T$, e.g., $T=60$, results in high SNR, but the temporal resolution is compromised as the ergodic assumption may be violated. For example, $T=15$ still works well for the single-scattering component as shown in the first panel in Fig. 3 (a)), in contrast to the significant degradation in the multiple-scattering contrast image, as shown in the first panel in Fig. 3 (h)). $T=30$ was thus used in the following functional CBF imaging as a trade-off between the SNR and the temporal resolution.

4.4 Steady-state CBF imaging and validation by fluorescence imaging

Figure 4 (c $\sim$ e) shows the contrast images corresponding to the hybrid, single- and multiple- scattering components. With the recommended values of $N=9$ and $Q=30/9$, the multiple-scattering contrast image showed the highest SNR, which makes the deep vasculature more visible than in the single-scattering or hybrid images (see white arrow in Fig. 4 (f)). The existence of this vessel was confirmed by the fluorescence image (the first image in Fig. 4 (f)) and the U-shaped contrast values crossing the vessel (Fig. 4 (g)).

 

Fig. 4. RMT-based separation in an in vivo CBF imaging experiment. (c $\sim$ e) corresponding to hybrid, single- and multiple-scattering contrast images, respectively; (a) the enlarged blue circle areas in (c $\sim$ e) and the corresponding fluorescence image showing the superficial vasculature details in (d). (b) and (g) show the contrast value curves crossing the selected vessel, indicated by the white arrows in (a) and (f). (f) The enlarged red circle areas in (c $\sim$ e) show more deep tissue vasculature details in (e) than in (c, d). The corresponding fluorescence image confirmed the presence of this vessel.

Download Full Size | PDF

The single-scattering contrast image also showed more vasculature details in the superficial layers (thicknesses: $<l_t/2$). A small vessel in Fig. 4 (a) (see the white arrow) was only visible in the single-scattering contrast image, not in the hybrid or multiple-scattering contrast images, or even in the fluorescence images. The existence of this vessel was confirmed by the U-shaped contrast values crossing the vessel (Fig. 4 (b)). The visibility of the small vasculature in the single-scattering contrast image was higher than in the fluorescence image.

4.5 Functional CBF responses in the superficial and deep layers

Figure 5 shows the CBF responses in the superficial and deep layers induced by hind paw stimulation. Figure 5 (a) and (e) show the single- and multiple-scattering contrast images at their maximum responses to hind paw stimulation with the optimal settings ($S=3$ and $T=30$), respectively. In addition to the improved SNR, the multiple-scattering contrast image was immune to specular noise (see the white arrow in Fig. 5 (a)). To compare the different responses in the superficial and deep layers of the rat cerebral cortex, we show the typical blood flow responses in four rectangular areas covering both the vasculature (blue boxes, Area1 and Area2) and tissue (red boxes, Area3 and Area4), as shown in Fig. 5 (a) and (e). The two separated tissue areas (Area3 and Area4) exhibit similar responses, and the responses in the deep layers (manifested by multiple-scattering components) are always larger than those in the superficial layer (manifested by the single-scattering components). To further compare the responses of the superficial and deep tissue layers, we segmented the major vessels using an automatic transfer learning-based method and divided the tissue area into grids of $7 \times 7$ pixel windows. A window with CBF changes in the contrast values of more than (10%) was marked as a significant location. Figure 5 (d) shows the significant locations in the deep (red points) and superficial layers (blue), respectively, with an odds ratio of 100:44, showing more significant hindpaw stimulation-induced brain activities in the deep layers. Previous studies have also confirmed that the coupling between functional stimulation-induced neural activities and tissue perfusion was primarily in the deep layer of the cortex [36,37].

 

Fig. 5. The proposed strategy is applied to the functional responses of blood flow and tissue perfusion to the electrical hind paw stimulation. (a) and (e) are single- and multiple-scattering contrast images, respectively. Four areas of interest covering the vasculature and tissue are indicated by the red and black boxes in (a) and (e). The corresponding responses of the blood flow and tissue perfusion estimated from the single- and multiple-scattering components are shown in (b), (c), (f) and (g), respectively. The responses by the hybrid contrast data are also plotted. (d) shows the locations in the tissue area with a significant change (>10%) based on the average of the baseline level in both the superficial (blue markers) and deep layers (red markers). (h) shows the locations in the vasculature area with a significant difference (>10%) in the average responses.

Download Full Size | PDF

Within both ROIs in Area1 and Area2, the vessel diameters were greater than $l_t/2$, implying that they crossed both the superficial and deep layers. Figure 5 (f) and (g) shows distinct responses in the single- and multiple-scattering components in both vessels. Vessels with more volume in the superficial layer had greater responses in the single-scattering component. Using a similar process as for Fig. 5(d), we compared the differences between the vessels in the superficial and deep layers, as shown in Fig. 5 (h), showing an odds ratio of 100:65 for >10% CBF changes in the deep layer (red points). Some overlapping red and blue points were observed in Fig. 5(h), which might reflect non-Newtonian characteristics due to the aggregated motion of red blood cells in in vivo CBF [38,39].

5. Discussions

Traditional LSCI is a mixture of both single- and multiple-scattering components; as a result, the deep and superficial layers interfere in the vasculature. The separation of the blood flow signals in the superficial and deep layers of the sample provides a unique opportunity to investigate stimulation-induced blood flow and perfusion changes, which is particularly important for tissues with layered structures, e.g., the cerebral cortex and the retinal fundus. The responses in different layers may reveal more specific physiological and pathological information.

When applied to other imaging modalities that use full-field coherent illumination, our method also eliminates interference from the superficial layer, allowing for more accurate measurements. For example, our skin has high scattering and thus contributes to the majority of the backscattered lights in the recorded signals, which biases the measurement of the deeper layers. A pair of polarizers can be used to extract the cross-polarization part since most superficial components have parallel polarization characteristics. In most cases, high illumination power is required to record sufficient signals, which is detrimental to the skin. Polarizers, however, could indiscriminately filter many trajectories from the deeper layers, resulting in insufficient deep tissue sampling. Furthermore, in many clinical applications, e.g., endoscopy and surgical microscopy, polarizers are always the second choice due to the limited work space.

In this study, a Wishart random matrix was constructed based on independently captured dynamic speckle images, with the imaging speed limited by temporal sampling. To improve the imaging speed, we can implement various strategies, such as multiple cameras and/or multiple illumination directions. Multiple-camera devices have been used in a variety of computer vision applications. Because speckle images from different cameras are treated as independent samples of the same ensemble, the data collection procedure can be significantly accelerated when multiple illumination directions are used. Nevertheless, similar verification and optimization procedures, such as EMC simulations and in vitro simulations, should be implemented to ensure that the separation of single- and multiple-scattering components is correct. The current Monte Carlo simulation applied $l_t=2.066 \mu m$ to minimize the computer memory consumption in thick tissue simulation while satisfying the Nyquist sampling theorem. Although $l_t$ in current Monte Carlo simulation is different from that of brain tissue or Intralipid phantom, the validity of current theoretical framework does not change.

Although the multiple-scattering contrast image revealed more vasculatures, the imaging depth was not substantially improved. The first reason for this is that out-of-phase suppression occurs in longer trajectories. The suppression is exponentially related to the path lengths of the multiply backscattered lights, limiting the fundamental penetration depth of LSCI. The second reason is absorption in biological tissue, which completely eliminates some light paths, especially in long trajectories, further limiting the imaging depth. Furthermore, the multiple-scattering contrast image exhibits blurring effects. More sophisticated deblurring methods are needed to further improve blood flow imaging in deep layers.

Although the speckle images could be successfully separated into single- and multiple- scattering components corresponding to layers at different depths, the exact depth information remained unresolved. This is an intrinsic limitation of LSCI as a $2 \mathrm {D}$ full-field imaging modality. However, our theoretical framework can be applied to other $3 \mathrm {D}$ imaging methods, e.g., multifocal reconstruction, DOT/DCT and OCT, where the separation of the single- and multiple-scattering components can be resolved with accurate depth coordinates. Other RM construction strategies can also be developed based on specific illumination and detection principles.

6. Conclusion

In conclusion, we investigated the spectral density of a Wishart random matrix constructed from dynamic speckle images and proposed a method for separating the single-scattering and multiple-scattering components. We used electric field Monte Carlo (EMC) simulations and an in vitro phantom experiments to verify the theory and optimize the model parameters. The proposed strategy was used to image blood flow in vivo. The separated single-scattering contrast image mainly showed superficial vasculature details, while the multiple-scattering contrast image revealed more information in deep tissue layers. We also demonstrated a distinct functional response under electrical hind paw stimulation in the superficial and deep layers of rat’s cerebral cortex. Our method extended the imaging ability of traditional LSCI in a wide range of applications.