1. Introduction

Structural and functional brain imaging has become a very important tool in neuroscience research. As a simple two-dimensional imaging method, laser speckle imaging (LSI) [1] has been widely used in imaging the cerebral vasculature and cerebral blood flow (CBF) [2], e.g. the vasculature in brain tumor [3] and retina [4], the functional changes of CBF under focal cerebral ischemia [2], hypotension [5] and/or peripheral electrical stimulation [6]. Both spatial [7] and temporal [8] laser speckle contrast analysis methods, i.e. spatial LASCA and temporal LASCA respectively, have been used to estimate the contrast values ($K^2$) which is inversely proportional to the relative velocity (v) [9]. However, spatial LASCA decreases the spatial resolution while temporal LASCA is under the cost of temporal resolution. Furthermore, the inevitable noise in both LASCA methods not only decreases the distinguishability of vessels, but also results in errors in estimating functional changes of blood flow.

In this paper, we develop a random process theory for explaining the speckle phenomena and derive the relation between the velocity of CBF and the probability distribution of speckle’s random process. With the help of the random process theory, both LASCA methods are directly derived as the spatial and temporal estimators, and the noise is also modeled and analyzed. Finally, based on the property of noise, a new random process estimator is proposed to improve the signal-noise ratio (SNR) of the estimation, which thus leads to a better structural and functional imaging of CBF.

2. Theory

2.1 Background and preparation

When coherent laser light illuminates a surface of rat’s cortex, the random interference patterns, i.e. laser speckles, are produced by the coherent addition of scattered laser lights with slightly difference in light-path lengths. The motion of scattering particles, i.e. the red cells in blood vessels, changes the temporal pattern of speckles. The principle of LASCA is to estimate the relative velocity of CBF by analyzing the changing speckle patterns.

Essentially, light is a kind of electromagnetic radiation. According to the electromagnetic theory, the area under coherent laser illumination can be modeled as a ‘complex electric field’ $E(t),t\in [0,+\infty )$ [10,11]. Following the dynamic light scattering (DLS) theory [12], the scatters’ velocity$v(t)$is related to$E(t)$’s autocorrelation function$g_1(\tau )$:

2.2 Random process theory for laser speckle phenomenon

Due to the randomness property of speckle phenomena, the intensity$I(t),t\in [0,+\infty )$of one speckle is considered as a single realization of the continuous intensity random process$X(t)$. In this section, we are going to develop the random process theory to explain the laser speckle phenomenon and derive the relation between$X(t)$and velocity$v(t)$.

At any time$t_1$, $I(t_1)$is a single realization of the corresponding intensity random variable $X(t_1)$ whose probability distribution is related to the velocity$v(t_1)$. In practice, a CCD camera is used to acquire the intensity information of the imaging plane. By carefully adjusting the f number of the imaging system, we can make the pixel size of the recorded image equal to the speckle size so that there is only one speckle in one pixel. Therefore, the ‘intensity’ value of any pixel is actually equal to the integrated intensity over the exposure durationT:

Suppose the exposure time T is designed to be very short ($~5ms$), then at any time$t_1$, the velocity $v(t^{\prime }),t^{\prime }\in [t_1,t_1+T]$ can be approximately considered as a constant. Thus, the distributions of random variables $\{X(t^{\prime }),t^{\prime }\in [t_1,t_1+T]\}$ are the same, and each time-limited continuous random process $X(t^{\prime }),t^{\prime }\in [t_1,t_1+T]$ can be considered as a wide-sense stationarity (WSS) random process. The properties of WSS random process [16] include: (a) its mean value ${\mu }_X(t^{\prime })$ is a constant when$t^{\prime }\in [t_1,t_1+T]$; (b) its temporal average value of the product of any two samples, i.e. $X(t^{\prime })$and$X(t^{″})$, $t^{\prime },t^{″}\in [t_1,t_1+T]$, depends only upon the time interval $\tau =t^{″}-t^{\prime }$.

Based on the property (a), when the number of ensemble samples$N\to +\infty$, we have:

According to the Siegert relation [12], we have:

Suppose we obtain infinite temporal realizations of$\{X(t^{\prime }),t^{\prime }\in [t_1,t_1+T]\}$, then we have:

By simply manipulating Eq. (6), we conclude:

Using$li{m}_{N\to +\infty }{〈Y(t_1)〉}_{ensemble}\to {\mu }_Y(t_1)$and$li{m}_{N\to +\infty }{〈Y^2(t_1)〉}_{ensemble}-{〈Y(t_1)〉}_{ensemble}^2\to {\sigma }_Y^2(t_1)$, we rewrite Eq. (7) and define the contrast of integrated intensity random variable$Y(t_1)$, i.e. $K_Y^2(t_1)$, as follows:

Based on different models describing the relation between $v(t)$ and $g_1(\tau )$, the relation between $K_Y^2(t)$ and $v(t)$ can also be derived directly [13,15,19,20]. For any single speckle, if we can obtain its ensemble realizations (samples) $\{{\tilde{I}}_i^2(t_1)\}$of the integrated intensity random variable $Y(t_1)$ at time $t_1$, the ensemble samples can be used to estimate the $K_Y^2$ and then the relative velocity v can be estimated.

2.3 The spatial and temporal estimators and the estimation noise

In practice, at any time$t_1$, we can only get one realization of $Y(t_1)$, i.e. $\tilde{I}(t_1)$, for each speckle (pixel). So we have to estimate $K_Y^2(t_1)$ in other ways. An alternative method is to estimate the ${\mu }_Y(t_1)$ and ${\sigma }_Y(t_1)$ using the intensities of the pixel ($I_0(t_1)$) and its surrounding pixels ($I_j(t_1),j=1,\dots ,M-1$) based on the assumption that velocities at the surrounding pixels are very close. When this assumption is true, the distributions of integrated intensity variable $\{Y_j(t_1)\}$ at different pixels are the same. Thus, spatial samples $\{{\tilde{I}}_j(t_1),j=0\dots M-1\}$ can be considered as ensemble samples approximately: ${〈Y(t_1)〉}_{spatial}=\frac{1}{M}{\displaystyle {\sum }_{j=0}^{M-1}{\tilde{I}}_j{t}_1}\approx {〈Y(t_1)〉}_{ensemble}$ and ${〈Y^2(t_1)〉}_{spatial}=\frac{1}{M}{\displaystyle {\sum }_{j=0}^{M-1}{\tilde{I}}_j^2(t_1)}\approx {〈Y^2(t_1)〉}_{ensemble}$. The estimations of ${\mu }_Y(t_1)$, ${\sigma }_Y(t_1)$ and $K_Y^2(t_1)$ can be obtained as follows:

It is noted that this Eq. (9) is actually the spatial LASCA method [7]. Such estimation works well in most in-vitro simulation studies when the moving scattering particles in a large area have the same velocity. However, this assumption does not hold for a complex biological tissue, such as brain, where there is: 1) an abundance of scattering particles such as red cells, and 2) a large number of vessels with different diameters and velocities.

Next, we analyze the estimation noise in spatial LASCA. Generally, the estimation of ${\mu }_s(t_1)$ and ${\sigma }_s(t_1)$ in a spatial window centered at the interesting pixel could be affected by (i) limited size of the window, which results in the estimation noises $N{A}_{{\mu }_s}(t_1)$ and$N{A}_{{\sigma }_s}(t_1)$; and (ii) the spatial inhomogeneity of velocities within the window, which could result in estimation noises $N{B}_{{\mu }_s}(t_1)$ and $N{B}_{{\sigma }_s}(t_1)$. Therefore, ${\mu }_s(t_1)$ and ${\sigma }_s(t_1)$can be modeled as follows:

Although a larger window (e.g. $5\times 5$or larger) can result in smaller$N{A}_{{\mu }_s}$and$N{A}_{{\sigma }_s}$based on central limit theory, the $N{B}_{{\mu }_s}$and$N{B}_{{\sigma }_s}$will increase due to higher inhomogeneity of velocities. In contrast, a smaller window (e.g. $3\times 3$) will result in smaller$N{B}_{{\mu }_s}$and$N{B}_{{\sigma }_s}$, but larger$N{A}_{{\mu }_s}$and$N{A}_{{\sigma }_s}$.

Another problem of spatial LASCA is the loss of spatial resolution. When we need a high spatial resolution image of cerebral vessels, spatial LASCA is not appropriate. To preserve the spatial resolution, the temporal estimator was developed.

If we assume that the velocity in the interesting pixel is constant during the time interval$[t_1,t_2]$, then the probability distributions of integrated intensity variables $\{Y(t^{\prime }),t^{\prime }\in [t_1,t_2]\}$ are the same. In practice, we use a CCD camera to record the temporal samples ($\{\tilde{I}(t^{\prime })\}$) of integrated intensity random process$Y(t^{\prime })$:

Again the temporal estimator $K_{temporal}^2$ is actually the temporal LASCA method [8]. The temporal LASCA is also contaminated by the estimation noise, i.e. the noise components $NA$ and $NB$ in ${\mu }_t(t_1)$ and ${\sigma }_t(t_1)$ as follows:

Temporal LASCA preserves the spatial resolution because ${\mu }_t(t_1)$ and ${\sigma }_t(t_1)$ are calculated based on $\tilde{I}(t_1+n\Delta t)(n=0,\dots ,N-1)$ of the same pixel. Therefore, high resolution structural imaging of cerebral blood vessel network is achieved. When imaging rat’s CBF under stable condition, the velocity in each speckle is approximately stable during the recording, and the noise component $NB$ can be ignored. However, if the blood velocity is changing, e.g. elicited in response to functional stimulation, the $NB$ noise will be prominent.

To achieve a high performance of statistical calculation (low $NA$ noise), usually more than 50 frames are used in the temporal LASCA. In practice, the frame rate of the CCD camera is limited, for example $10fps$ in this study. So, the temporal resolution is compromised. Fewer frames, i.e. small time window, will reduce the loss of temporal resolution, but leads to larger $NA$ noise.

2.4 Random process estimator

Based on the properties of integrated intensity random process$Y(t)$, we propose a new estimator, called random process estimator hereafter, which provides high SNR estimation and high spatiotemporal resolution for both structural and functional imaging.

The velocity of CBF can be considered as a constant within a short exposure duration (e.g. $~5ms$). For an interesting pixel, we calculate the $\{{\sigma }_s(t_1+n\Delta t),n=0,\dots ,N-1\}$ using a $3\times 3$ window from each frame. Because the window is small, there would be large noise component $N{A}_{{\sigma }_s}$ and small $N{B}_{{\sigma }_s}$ in $\{{\sigma }_s(t_1+n\Delta t),n=0,\dots ,N-1\}$. Then, Eq. (11) is simplified into: ${\sigma }_s(t_1+n\Delta t)\approx {\sigma }_Y(t_1+n\Delta t)+N{A}_{{\sigma }_s}(t_1+n\Delta t)$, where$n=0,\dots ,N-1$.

It is noted that each ${\sigma }_Y(t_1+n\Delta t)$ is related to the velocity $v(t_1+n\Delta t)$ and thus could change during the imaging procedure. The noise component $N{A}_{{\sigma }_s}$ is due to the limited number of samples, so it is affected by the number of samples. Based on the central limit theorem, $N{A}_{{\sigma }_s}$can be hypothesized to follow a zero-mean Gaussian distribution, i.e. $N(0,{\sigma }_1^2)$, which will be confirmed based on the simulation data in Discussion. Therefore, the discrete sequences $\{N{A}_{{\sigma }_s}(t_1+n\Delta t),n=0,\dots ,N-1\}$ can be considered as independent and identically-distributed (IID) random variables with the same Gaussian distribution.

Then we calculate the average of $\{{\sigma }_s(t_1+n\Delta t)\}$ as the random process estimator for averaged${\overline{\sigma }}_Y(t_1)$:

${\mu }_t(t_1)$is used to estimate the averaged ${\overline{\mu }}_Y(t_1)$ in the random process estimator, so there is no loss of spatial resolution in the calculation of${\mu }_{rpe}(t_1)$. Actually, each $\tilde{I}(t_1+n\Delta t)$ can be considered as a ‘spatial estimation’ of ${\mu }_s(t_1+n\Delta t)$ with a $1\times 1$ spatial window, which leads to a large noise component $N{A}_{{\mu }_s}$ without $N{B}_{{\mu }_s}$ noise. Ignoring the $NB$ noise in Eq. (10), we have: $\tilde{I}(t_1+n\Delta t)={\mu }_Y(t_1+n\Delta t)+N{A}_{{\mu }_s}(t_1+n\Delta t)$. Then we get:

Similarly, we model $N{A}_{{\mu }_s}(t_1+n\Delta t)$ as a sequence of IID random variables with a zero-mean Gaussian distribution based on the central limit theorem, i.e. $N(0,{\sigma }_2^2)$ (this assumption will also be confirmed based on the simulation data in Discussion). Therefore, when $N\to +\infty$, $\frac{1}{N}{\displaystyle {\sum }_{n=0}^{N-1}N{A}_{{\mu }_s}(t_1+n\Delta t)}\to 0$ and then${\mu }_{rpe}(t_1)={\overline{\mu }}_Y(t_1)=\frac{1}{N}{\displaystyle {\sum }_{n=0}^{N-1}{\mu }_Y(t_1+n\Delta t)}$, where ${\overline{\mu }}_Y(t_1)$ corresponds to the averaged velocity $\overline{v}$ within$[t_1,t_1+(N-1)\Delta t]$.

Finally, the new random process estimator is used to estimate $K_Y^2(t_1)$ for the averaged velocity $\overline{v}$ during time interval $[t_1,t_1+(N-1)\Delta t]$ as follows:

Compared with spatial LASCA, $K_{rpe}^2$provides much higher spatial resolution because the calculation of ${\sigma }_{rpe}$ uses only a $3\times 3$ window. Furthermore, because the frame number N only affects the averaging (denoising) performance of the random process estimator, we can use fewer frames (e.g. $N=10$) than the temporal LASCA ($N=50$images) to achieve a higher temporal resolution. Even using only 10 frames, the denoising performance of random process estimator is still significant while temporal LASCA is contaminated by both $NA$ and$NB$. Therefore, random process estimator is very appropriate for functional imaging where both spatial and temporal resolutions need to be preserved as much as possible.

Although temporal LASCA can provide high resolution image of blood vessel network based on a large number of raw images, the noise component $NA$ is unavoidable. Furthermore, some physiologic changes, e.g. body temperature, heart rate, blood pressure and medicine, could also lead to a time-varying CBF and thus the $NB$ noise, which will decrease the distinguishability of small vessels. On the other hand, larger N in random process estimator will have better denoising performance according to Eq. (16) and Eq. (17). Therefore, although random process estimator produces a slight loss in spatial resolution compared to temporal LASCA, it significantly improves the distinguishability of vessels.

3. Experiments

3.1 Imaging setup

The imaging plane was illuminated with a $632nm$ He-Ne laser beam source (1508P-O, Uniphase, CA) which was reshaped by optical lens to expand the range of illumination. A 12-bit cooled CCD camera (Sensicam SVGA, Cooke, MI) with a $60mm$ macro (1:1 maximum reproduction ratio) $f/2.8$lens was held and focused on the imaging plane. Exposure time of the CCD camera was set to $5ms$. Frame rate was $10fps$. The resolution of image recorded by the CCD camera was $704\times 704$ pixels corresponding to an imaging area about $4.7\times 4.7m{m}^2$.

3.2 In-vitro simulation experiment

A phantom is created to mimic blood vessels in vitro by flowing blood through plastic tubing using a syringe pump (PHD2000, Harvard Apparatus, MA). Real blood was obtained by ex-sanguination of rats destined for sacrifice and mixed with heparin sodium (10 IU/ml, Sigma Chemical Co., St. Louis, MO) to prevent clotting. A syringe filled with such blood was fitted into the syringe pump and connected to a polyethylene tube (internal diameter$0.034"$, PE90). Blood was infused into the tube at $2mm/s$. When the flow was steady after $~5min$, the tube was imaged using the imaging setup described above. 600 frames were acquired for analysis.

3.3 In-vivo structural and functional LSI experiment

The experimental protocol used in this study has been approved by the Animal Care and Use Committee of the Johns Hopkins Medical Institutions. The female Sprague-Dawley rats ($~320g$) were used in this study. To avoid the influence of surgical preparation on the rats’ CBF, the thinned skull preparation was done on the first day. On the second day, the rats were taken back for the structural and functional LSI experiment.

On the first day, the rat was anesthetized with intraperitoneal (IP) injection of mixture of $90mg/kg$ of Ketamine and $10mg/kg$ of Xylazine. The animal was placed in a stereotactic frame (David Kopf Instruments, Tujunga, CA) throughout the experiment. The scalp was shaved and disinfected with $\text{7}0\%$ ethanol and povidone-iodine solution. All procedures were performed using standard sterile precautions. After a midline scalp incision, the galea and periosteum overlying the parietal bone bilaterally were swept and retracted laterally. A $5mm\times 5mm$ area centered at $2.5mm$ lateral and $2.5mm$ posterior to the bregma over the right cortex was thinned using a high speed dental drill (Fine Science Tools Inc. North Vancouver, Canada), until the inner cortical layer of bone was encountered. Rectal temperature was maintained at ${37}^{°}\text{C}$ during the experiment using a homeothermic blanket system (model TP-500 T/Pump, Gaymar Industries, Inc., USA).

On the second day, the rat was first held in a transparent chamber with $3\%$ isofluorane gas and room air flow until becoming drowsiness. Its mouth and nose were then placed within a customized anesthesia mask in a well-fitting rodent size, which was connected to a C-Pram circuit designed to deliver and evacuate the gas through one tube. A mixed flow of $1.5\%$ isofluorane, $80\%$oxygen and room air was delivered to the mask at the flow of $2L/min$ for anesthesia. 600 frames were recorded for the structural imaging experiment.

Afterwards, subcutaneous needle electrode pairs (Safelead F-E3-48, Grass-Telefactor, West Warwick, RI) were inserted just under the skin of left hindpaw, one on the top of the hindpaw between 2 and 3 digits, the other one on the back skin of the hind paw, without direct contact to the nerve bundle. An isolated constant current stimulator (DS3, Digitimer Ltd., Hertfordshire, England) was used for the electrical stimulation of the hindpaw. Positive current pulses of $3.5mA$ magnitude and $200\mu s$ duration at a frequency of $3Hz$ were used for hind paw stimulation. Each stimulation trial consisted of pre-stimulus 200 frames ($20s$), 200 stimulation frames ($20s$), and 400 post-stimulus frames ($40s$), which was repeated ten times with several-minute in-between breaks.

4. Results

4.1 The denoising performance of random process estimator

The denoising performance of random process estimator is evaluated in the in-vitro simulation experiment. During the experiment, the blood flow velocity in the polyethylene tube was strictly maintained at $2mm/s$ using the syringe pump. According to Eq. (16) and Eq. (17), the noise parts $\frac{1}{N}{\displaystyle {\sum }_{n=0}^{N-1}N{A}_{{\mu }_s}(t_1+n\Delta t)}$ and $\frac{1}{N}{\displaystyle {\sum }_{n=0}^{N-1}N{A}_{{\sigma }_s}(t_1+n\Delta t)}$ in the random process estimator would be very closed to zeros with$N=600$. Therefore, for any pixel, ${\mu }_{rpe}$and ${\sigma }_{rpe}$ based on all 600 raw images were considered as the true ${\mu }_Y$ and${\sigma }_Y$. Figure 1(a) shows one frame of the 600 raw images. Figure 1(c) is the estimated contrast image $K_Y^2$ using random process estimator. One pixel (the white point in the white box area in Fig. 1(c)) is selected arbitrarily to compare the estimation results based on spatial LASCA, temporal LASCA and random process estimator.

 

Fig. 1 The denoising performance of random process estimator: (a) one frame of raw images; (b) the estimation results of ${\sigma }_Y$ using different methods; (c) the contrast image $K_Y^2$ estimated based on all 600 raw images using the random process estimator; (d) the estimated $K_Y^2$ using different methods.

Download Full Size | PDF

For the selected pixel in Fig. 1(c), the result of spatial LASCA is calculated using a $7\times 7$ spatial window at each second ($t_1=0,1,\dots ,59s$) while the results of temporal LASCA and random process estimator are calculated using a $1s$ temporal window (10 frames) starting from each second. Figure 1(b) shows the estimations of ${\sigma }_Y^2$ based on spatial LASCA, temporal LASCA and random process estimator respectively. As expected, all estimations vary around the true value (white dashed line). Among all three methods, temporal LASCA produces the largest fluctuations because of large noise component $NA$ (small temporal window) in Eq. (15). Although both LASCA methods are contaminated by $NA$ noise, spatial

LASCA still provides a better estimation using a comparably bigger window size (spatial LASCA: $N=7\times 7=49$; temporal LASCA: $N=10$). Furthermore, according to the hydrodynamics, the blood velocity is still not distributed uniformly within the tube. Therefore, besides the $NA$ noise, there is also $NB$ noise in the result of both LASCA methods. Compared to the conventional LASCA methods, random process estimator provides the best estimation of ${\sigma }_Y$ in the sense of smallest fluctuations. Since the estimations of ${\mu }_Y$ are the same in temporal LASCA and random process estimator, we did not show the results of${\mu }_Y$.

Figure 1(d) shows the estimations of $K_Y^2$ at each second ($t_1=0,1,\dots ,59s$) based on spatial LASCA, temporal LASCA and random process estimator where random process estimator outperformed the others again. To quantitatively compare the denoising performances of different methods, we define the signal-noise ratio (SNR) for the estimation of $K_Y^2$ as follows:

Although spatial LASCA provides a slightly better estimation than temporal LASCA, the loss of spatial resolution is unavoidable and usually unacceptable. Since small vessels play an important role in the in-vivo structural and functional studies, we’ll only compare temporal LASCA and random process estimator in the rest sections.

4.2 Random process estimator for structural imaging

In the in-vivo experiment, 600 raw images ($60s$) are recorded for imaging the structure of blood vessel network. During the imaging procedure, the rats were anesthetized and kept in stable condition, so blood velocities in the imaging area are assumed to be stable. Figure 2 shows the contrast images estimated from the first 2 (Fig. 2(a,d)), 10 (Fig. 2(b,e)) and 80 (Fig. 2(c,f)) frames using random process estimator and temporal LASCA respectively. For either method, the distinguishability of small vessels is improved when more frames are used. For example, the small vessel designated with black arrow in Fig. 2 (b) and (c) is hard to see in Fig. 2(a).

 

Fig. 2 Structural imaging of rat’s cerebral blood vessels: (a, b, c) show the contrast images $K_Y^2$ estimated from the first 2, 10 and 80 frames using random process estimator; (d, e, f) show the contrast images $K_Y^2$ estimated from the first 2, 10 and 80 frames using temporal LASCA.

Download Full Size | PDF

According to Eq. (15), when more frames are used, the noise component $NA$ decreases but still exists in the result of temporal LASCA. Furthermore, if there are small changes of blood velocity during the imaging procedure, the $NB$ noise is also introduced in. However, using random process estimator, the noise part is rapidly averaged out according to Eq. (16), 17) provided that more frames are used. More precisely, based on the same number of frames, random process estimator always provides better estimation and thus better distinguishability than temporal LASCA. For example, the small vessels in the white circled area in Fig. 2(c) are more distinguishable than in Fig. 2(f).

The contrast values along the vertical white line in Fig. 2(c) estimated from different number of frames using the two methods is plotted in Fig. 3 . There are six vessels denoted as V1~V6 crossing the line. In Fig. 3, the results of random process estimator are always smoother than temporal LASCA. In particular, based on only 2 frames, the noise in the result of temporal LASCA was prominent (Fig. 3(a)), while the result of random process estimator still reveals the denoising aspect especially in the vessel areas V1 and V2. For 10 frames (Fig. 3(b)), the random process estimator already provides an efficient estimation of $K_Y^2$ with high

 

Fig. 3 Denoising performance of random process estimator for structural imaging: the contrast values $K_Y^2$ along the vertical line in Fig. 2 (c) estimated from different number of frames: N=2 (a), N=10 (b), N=80 (c) using random process estimator and temporal LASCA.

Download Full Size | PDF

SNR while the result of temporal LASCA is still very noisy. Furthermore, the denoising performance of random process estimator is better in the tissue areas compared to vessel areas.

Although more frames can be used in the calculation by random process estimator, no significant improvement is found compared to Fig. 2(c). Therefore, 80 images are enough for the structural imaging using random process estimator. In conclusion random process estimator provides better distinguishability of structural information of cerebral blood vessels.

4.3 Random process estimator for functional imaging

In the functional stimulation experiment, each trial includes recording 800 raw images ($\text{8}0s$) continuously, i.e. 200 pre-stimulus frames ($20s$), 200 stimulation frames ($20s$) and 400 post-stimulus frames ($40s$). Images within each second (10 frames) are used to calculate one contrast image using random process estimator and temporal LASCA. Therefore, for each trial, we obtain 80 contrast images ($K_{-20}^2,\dots ,K_{59}^2$).$K_{-20}^2$is used to normalize the change of contrast values as follows ($n=-19,\dots ,59$):

Because the stimulation procedure is 20s, other cortex areas may also be activated besides the somatosensory cortex. In this study, the area covering the activated somatosensory cortex is selected as the region of interest in the following discussions.

To analyze the denoising performance for functional imaging, the $\{K_{-20}^2,F_{-19},\dots ,F_{-16}\}$ of the interesting region in pre-stimulation stage of the first trial are shown in Fig. 4 using random process estimator (the first row) and temporal LASCA (the second row) respectively. Since the rat’s condition was kept stable before the stimulation, there should not be significant changes in the contrast images and thus the values in images $\{F_n,n=-19,\dots ,-1\}$ are expected to be much closed to zero. However, as shown in Fig. 4, temporal LASCA leads more noises than random process estimator. To quantify this estimation error, we define the ‘averaged errors’ in each $F_n$ as follows:

 

Fig. 4 Spatiotemporal change of $K_Y^2$ during the pre-stimulation stage: the first row shows the $K_{-20}^2,F_{-19},\dots ,F_{-16}$ (from the left to the right) using random process estimator; the second raw shows the $K_{-20}^2,F_{-19},\dots ,F_{-16}$ (from the left to the right) using temporal LASCA.

Download Full Size | PDF

Figure 5 shows the averaged errors of single trial estimation using temporal LASCA and random process estimator respectively. For random process estimator, the averaged errors in all $\{F_n,n=-19,\dots ,-1\}$ is only $0.31\pm 0.03$ ($mean\pm s.t.d.$), in contrast with $1.36\pm 0.09$ for temporal LASCA.

 

Fig. 5 The averaged errors of single trial result using temporal LASCA and random process estimator respectively.

Download Full Size | PDF

$\{F_0,\dots ,F_{59}\}$correspond to functional changes during the stimulation and post-stimulation procedure. By averaging the $\{F_0,\dots ,F_{59}\}$ of all ten trials, the averaged spatiotemporal changing map ($\{{\overline{F}}_0,\dots ,{\overline{F}}_{59}\}$) is obtained. Figure 6 and Fig. 7 show the averaged changing map$\{{\overline{F}}_0,\dots ,{\overline{F}}_{33}\}$, using random process estimator and temporal LASCA respectively. As a reference, the first frames in Fig. 6 and Fig. 7 are the structural image calculated by random process estimator.

 

Fig. 6 The changing patterns of CBF under functional stimulation calculated by random process estimator.

Download Full Size | PDF

 

Fig. 7 The changing patterns of CBF under functional stimulation calculated by temporal LASCA.

Download Full Size | PDF

In the functional imaging experiment, the blood flow in the somatosensory cortex changes in response to electrical stimulation. Therefore, besides the $NA$ noise, the $NB$ noise is also introduced into the result of temporal LASCA. Hence, the result of temporal LASCA (Fig. 7) is too noisy to precisely quantify the changes of CBF in small vessels. However, random process estimator (Fig. 6) provides high SNR and high spatiotemporal resolution for the functional changing pattern of CBF.

In Fig. 6, the response of CBF induced by functional stimulation starts from $4s$ and reaches maximum at$20s~22s$, then recovers to the baseline at$32s$. All significant changes can be located to small vessels (e.g. V1~V4 in Fig. 6). V1 is the vessel with the most significant changes, which is activated earlier than other vessels and shows longer recovery duration. With the help of high spatial resolution, we also find that not all vessels are activated in this area (e.g. V5).

One pixel in vessel V1 is chosen for quantitative analysis (the white circle in the first frame of Fig. 6). The averaged value and standard deviation of $\{F_n,n=-19,\dots ,59s\}$ at this pixel across all 10 trials are plotted in Fig. 8 for both methods. The result of temporal LASCA exhibits a noisy fluctuations ($-10\%~-40\%$) even in the pre-stimulation stage, so that only changes greater than $50\%$ can be considered as significant change during stimulation. However, the result of random process estimator shows much smaller fluctuations and a more stable profile. Furthermore, the standard deviation of the $\{F_n,n=-19,\dots ,59s\}$ over all 10 trials using random process estimator are always less than that using temporal LASCA (Fig. 8 (b)), implying that random process estimator also suppress the noisy variations over trials.

 

Fig. 8 The functional changes of contrast values corresponding to the pixel designated by a white circle in the first frame of Fig. 6: (a) the averaged changes, i.e. ${\overline{F}}_n,n=-19,\dots ,59s$, using temporal LASCA and random process estimator; (b) the standard deviation of the $F_n,n=-19,\dots ,59$s across all ten trials using temporal LASCA and random process estimator.

Download Full Size | PDF

5. Discussion

5.1 The Gaussian noise assumptions for random process estimator

In the development of random process estimator, the most important hypotheses are the Gaussian noise models in Eq. (16) and Eq. (17). Using the data acquired in the in-vitro experiment, we could test the hypotheses.

For the selected pixel in Fig. 1(c), 600 estimations of ${\sigma }_Y$ can be obtained using a $3\times 3$ spatial window. The noise component $N{A}_{{\sigma }_s}$ of each estimation is then calculated by subtracting the true value${\sigma }_Y$. The probability density function (pdf) is then calculated from all 600 $N{A}_{{\sigma }_s}$ by the kernel smoothing density estimation (Fig. 9(a) ). The pdf is very close to a Gaussian distribution with $0\pm 2.158\times {10}^{-3}$ ($mean\pm s.t.d.$). The P-P plot of all 600 $N{A}_{{\sigma }_s}$ in Fig. 9(b) is also very close to that of a Gaussian distribution. The hypothesis that the distribution of $N{A}_{{\sigma }_s}$ is a Gaussian distribution is further confirmed by the K-S test with a high z value ($z=0.859$).

 

Fig. 9 Gaussian assumptions confirmed by the data acquired in the in-vitro experiment: (a) the pdf of$N{A}_{{\sigma }_s}$; (b) the P-P plot of all 600 $N{A}_{{\sigma }_s}$s; (c) the pdf of$N{A}_{{\mu }_s}$; (d) the P-P plot of all 600 $N{A}_{{\mu }_s}$s.

Download Full Size | PDF

Similarly, the distribution of the noise component$N{A}_{{\mu }_s}$, i.e. ${\tilde{I}}_i-{\mu }_Y(i=1,\dots ,600)$, is also tested by K-S test with a high$z=0.848$. The pdf of the distribution ($0\pm 9.405\times {10}^{-3}$) and the P-P plot are shown in Fig. 9(c) and (d) respectively. These tests confirm that the distribution of $N{A}_{{\mu }_s}$ is also a Gaussian distribution.

5.2 Comparison with hybrid temporal and spatial method

Recently, a hybrid temporal and spatial method was proposed to improve the SNR by spatially averaging ($5\times 5$spatial window) the result of temporal LASCA which was successfully applied in imaging rat’s retinal blood flow [21]. Similar to spatial and temporal LASCA, the hybrid method works well when blood flow is constant within both spatial and temporal window, i.e. $NB$noise does not exist. We further process the same data set in our study using the algorithm in [21] to compare the performances of random process estimator and hybrid method in SNR, distinguishability of small vessels and revealing the functional responses in small vessels.

In the in-vitro simulation experiment, for the selected pixel in Fig. 1(c), the hybrid method produces a higher SNR (16.48) than temporal LASCA (2.75) and spatial LASCA (10.34), which, however, is still significantly lower than the SNR by random process estimator (31.09), because the use of a $5\times 5$spatial window increases the$NB$noise while suppressing the$NA$noise. In the structural imaging of the selected area in Fig. 10(a) , the details of small vessels (indicated with black arrows) would be hardly visible by the hybrid method with a$5\times 5$spatial window (Fig. 10(b)). Although using a $3\times 3$ spatial window in hybrid method slightly improves the spatial resolution (Fig. 10(c)), but the distinguishability of the vessels is still not as good as that by random process estimator (Fig. 10(d)) due to the decrease of SNR. Furthermore, the hybrid method cannot suppress the$NB$noise when there are significant blood flow changes within the temporal window. For example, $NB$noise is so large that the small functional changes are strongly contaminated by background noise in temporal LASCA (Fig. 11(b) ) under the functional stimulation. Hybrid method will filter out noise and the weak functional signals as well (Fig. 11(c)). However, random process estimator well reveals the functional changes in the small vessel while removing the noise (Fig. 11(d)).

 

Fig. 10 The vascular details in the white box area in (a) using hybrid method with$5\times 5$ (b) and$3\times 3$ (c) spatial window, random process estimator (d).

Download Full Size | PDF

 

Fig. 11 The functional CBF changes of selected area in (a) at the 8th sec estimated by temporal LASCA (b) as in Fig. 7, hybrid method (c) and random process estimator (d) as in Fig. 6.

Download Full Size | PDF

6. Conclusion

In this paper, we develop the random process theory for the laser speckle phenomenon and thus propose a new random process estimator to estimate the true value of $K_Y^2$. Compared with traditional spatial LASCA and temporal LASCA, random process estimator provides higher SNR and higher spatiotemporal resolution contrast images for both structural and functional imaging of cerebral blood vessels and flow.