Random process estimator for laser speckle imaging of cerebral blood flow

doi:10.1364/OE.18.000218

« Show journal navigation

Random process estimator for laser speckle imaging of cerebral blood flow

Peng Miao, Nan Li, Nitish V. Thakor, and Shanbao Tong »View Author Affiliations

Optics Express, Vol. 18, Issue 1, pp. 218-236 (2010)
http://dx.doi.org/10.1364/OE.18.000218

View Full Text Article

Acrobat PDF (1864 KB)

Journal Search
Article Lookup

Article Tools

Citations

Alert me when this article is cited
Export Citation/Save

Full Text
Article Info
References (21)
Cited By
Figures (11)
Metrics

Abstract

In this paper, we develop a random process theory to explain the laser speckle phenomena. The relation between the probability distribution of speckle’s integrated intensity random process $Y (t)$ and the relative velocity $v (t)$ is derived. Based on the random process theory, traditional spatial or temporal laser speckle contrast analysis (i.e. spatial or temporal LASCA) can be derived as the spatial or temporal estimators respectively. Both spatial LASCA and temporal LASCA suffer from noise due to insufficient statistics and nonstationarity in either spatial or temporal domain. Furthermore, either LASCA results in a reduction of spatial or temporal resolution. A new random process estimator is proposed and able to overcome these drawbacks. In an in-vitro study, random process estimator outperforms either spatial LASCA or temporal LASCA by providing much higher SNR (random process estimator vs. spatial LASCA vs. temporal LASCA: 33.64±6.87 ( $m e a n \pm s . t . d .$ ) vs. 9.08±2.85 vs. 3.83±1.05). In an in-vivo structural imaging study, random process estimator efficiently suppresses the noise in contrast image and thus improves the distinguishability of small vessels. In a functional imaging study of cerebral blood flow change in the somatosensory cortex induced by rat’s hind paw stimulation, random process estimator provides much lower estimation errors in single trial data (random process estimator vs. temporal LASCA: 0.31±0.03 vs. 1.36±0.09) and finally leads to higher resolution spatiotemporal patterns of cerebral blood flow.

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

J. Briers, “Laser Doppler, speckle and related techniques for blood perfusion mapping and imaging,” Physiol. Meas. 22(4), 35–66 (2001). [CrossRef]

] has been widely used in imaging the cerebral vasculature and cerebral blood flow (CBF) [2

A. K. Dunn, H. Bolay, M. A. Moskowitz, and D. A. Boas, “Dynamic Imaging of Cerebral Blood Flow Using Laser Speckle,” J. Cereb. Blood Flow Metab. 21(3), 195–201 (2001). [CrossRef] [PubMed]

], e.g. the vasculature in brain tumor [3

D. Zhu, W. Lu, Y. Weng, H. Cui, and Q. Luo, “Monitoring thermal-induced changes in tumor blood flow and microvessels with laser speckle contrast imaging,” Appl. Opt. 46(10), 1911–1917 (2007). [CrossRef] [PubMed]

] and retina [4

H. Cheng, Y. Yan, and T. Duong, “Temporal statistical analysis of laser speckle images and its application to retinal blood-flow imaging,” Opt. Express 16(14), 214–219 (2008). [CrossRef]

], the functional changes of CBF under focal cerebral ischemia [2

A. K. Dunn, H. Bolay, M. A. Moskowitz, and D. A. Boas, “Dynamic Imaging of Cerebral Blood Flow Using Laser Speckle,” J. Cereb. Blood Flow Metab. 21(3), 195–201 (2001). [CrossRef] [PubMed]

], hypotension [5

A. Kharlamov, B. R. Brown, K. A. Easley, and S. C. Jones, “Heterogeneous response of cerebral blood flow to hypotension demonstrated by laser speckle imaging flowmetry in rats,” Neurosci. Lett. 368(2), 151–156 (2004). [CrossRef] [PubMed]

] and/or peripheral electrical stimulation [6

T. Durduran, M. G. Burnett, G. Yu, C. Zhou, D. Furuya, A. G. Yodh, J. A. Detre, and J. H. Greenberg, “Spatiotemporal Quantification of Cerebral Blood Flow During Functional Activation in Rat Somatosensory Cortex Using Laser-Speckle Flowmetry,” J. Cereb. Blood Flow Metab. 24(5), 518–525 (2004). [CrossRef] [PubMed]

]. Both spatial [7

J. Briers and S. Webster, “Laser speckle contrast analysis (LASCA): a nonscanning, full-field technique for monitoring capillary blood flow,” J. Biomed. Opt. 1(2), 174–179 (1996). [CrossRef]

] and temporal [8

H. Cheng, Q. Luo, S. Zeng, S. Chen, J. Cen, and H. Gong, “Modified laser speckle imaging method with improved spatial resolution,” J. Biomed. Opt. 8(3), 559–564 (2003). [CrossRef] [PubMed]

] 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

P. Miao, M. Li, H. Fontenelle, A. Bezerianos, Y. Qiu, and S. Tong, “Imaging the Cerebral Blood Flow with Enhanced Laser Speckle Contrast Analysis (eLASCA) by Monotonic Point Transformation,” IEEE Trans. Biomed. Eng. 56(4), 1127–1133 (2009). [CrossRef] [PubMed]

]. 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

P. Lemieux and D. Durian, “Investigating non-Gaussian scattering processes by using nth-order intensity correlation functions,” J. Opt. Soc. Am. A 16(7), 1651–1664 (1999). [CrossRef]

,11

D. Boas and A. Yodh, “Spatially varying dynamical properties of turbid media probed with diffusing temporal light correlation,” J. Opt. Soc. Am. A 14(1), 192–215 (1997). [CrossRef]

]. Following the dynamic light scattering (DLS) theory [12

B. Berne, and R. Pecora, Dynamic Light Scattering: with Applications to Chemistry, Biology and Physics (Dover Publications, 2000).

], the scatters’ velocity

v (t)

is related to

E (t)

’s autocorrelation function

g_{1} (τ)

g_{1} (τ) = {〈 E (t) E^{*} (t - τ) 〉}_{t e m p o r a l} / {〈 E (t) E^{*} (t) 〉}_{t e m p o r a l} .

(1)

whereτis the delay time;

E^{*} (t)

is the complex conjugate of

E (t)

;

{〈 \cdot 〉}_{t e m p o r a l}

expresses the temporal average. The relation between

v (t)

and

g_{1} (τ)

relies on different models. For example, based on the assumption of Lorentzian distribution of velocity, Fercher and Briers [13

A. Fercher and J. Briers, “Flow visualization by means of single-exposure speckle photography,” Opt. Commun. 37(5), 326–330 (1981). [CrossRef]

] proposed the first model:

g_{1} (τ) = \exp (- τ / τ_{0})

(

τ_{0}

is the correlation time,

τ_{0}^{- 1} \propto v

); based on diffusing wave spectroscopy (DWS) model [14

D. J. Pine, D. A. Weitz, P. M. Chaikin, and E. Herbolzheimer, “Diffusing wave spectroscopy,” Phys. Rev. Lett. 60(12), 1134–1137 (1988). [CrossRef] [PubMed]

], Bandyopadhyay et al. [15

R. Bandyopadhyay, A. Gittings, S. Suh, P. Dixon, and D. Durian, “Speckle-visibility spectroscopy: A tool to study time-varying dynamics,” Rev. Sci. Instrum. 76(9), 093110 (2005). [CrossRef]

] suggested a more rigorous model:

g_{1} (τ) = \exp (- γ {(6 τ / τ_{0})}^{1 / 2})

(γ is a constant parameter,

τ_{0}^{- 1} \propto v

). In practice, it is very difficult to measure

E (t)

and

g_{1} (τ)

directly, but it is common to study the intensity

I (t), t \in [0, + \infty)

, defined as

I (t) = {| E (t) |}^{2}

. Generally, the recorded laser speckle image contains the intensity information of a composition of speckles with same speckle size

s = 2.44 λ f/ #

(λis the wavelength of laser light,

f/ #

is the f number of the imaging system) [2

A. K. Dunn, H. Bolay, M. A. Moskowitz, and D. A. Boas, “Dynamic Imaging of Cerebral Blood Flow Using Laser Speckle,” J. Cereb. Blood Flow Metab. 21(3), 195–201 (2001). [CrossRef] [PubMed]

]. For each speckle, the scattering particles’ velocity

v (t)

is related to

E (t)

and thus related to the statistical property of

I (t)

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:

\tilde{I} (t) = \frac{1}{T} \int_{t}^{t + T} I (t^{'}) d t^{'}

(2)

where

\tilde{I} (t)

is a single realization of integrated intensity random process

Y (t) = \frac{1}{T} \int_{t}^{t + T} X (t^{'}) d t^{'}

Suppose the exposure time T is designed to be very short (

~ 5 m s

), then at any time

t_{1}

, the velocity

v (t^{'}), t^{'} \in [t_{1}, t_{1} + T]

can be approximately considered as a constant. Thus, the distributions of random variables

{X (t^{'}), t^{'} \in [t_{1}, t_{1} + T]}

are the same, and each time-limited continuous random process

X (t^{'}), t^{'} \in [t_{1}, t_{1} + T]

can be considered as a wide-sense stationarity (WSS) random process. The properties of WSS random process [16

G. Grimmett, and D. Stirzaker, Probability and random processes (Oxford University Press, 2001).

] include: (a) its mean value

μ_{X} (t^{'})

is a constant when

t^{'} \in [t_{1}, t_{1} + T]

; (b) its temporal average value of the product of any two samples, i.e.

X (t^{'})

and

X (t^{″})

t^{'}, t^{″} \in [t_{1}, t_{1} + T]

, depends only upon the time interval

τ = t^{″} - t^{'}

Based on the property (a), when the number of ensemble samples

N \to + \infty

, we have:

μ_{Y} (t_{1}) \approx {〈 Y (t_{1}) 〉}_{e n s e m b l e} = {〈 \frac{1}{T} \int_{t_{1}}^{t_{1} + T} X (t^{'}) d t^{'} 〉}_{e n s e m b l e} = \frac{1}{T} \int_{t_{1}}^{t_{1} + T} {〈 X (t^{'}) 〉}_{e n s e m b l e} d t^{'} \approx μ_{X} (t_{1})

(3)

Based on the property (b), the ‘intensity autocorrelation function’

g_{2} (τ)

X (t^{'})

during the exposure

t^{'} \in [t_{1}, t_{1} + T]

exists:

g_{2} (τ) = {〈 X (t^{'}) X (t^{'} - τ) 〉}_{t e m p o r a l} / {〈 X (t^{'}) 〉}_{t e m p o r a l}^{2} .

(4)

According to the Siegert relation [12

B. Berne, and R. Pecora, Dynamic Light Scattering: with Applications to Chemistry, Biology and Physics (Dover Publications, 2000).

], we have:

{〈 X (t^{'}) X (t^{'} - τ) 〉}_{t e m p o r a l} = {〈 X (t_{1}) 〉}_{e n s e m b l e}^{2} (1 + β {| g_{1} (τ) |}^{2}),

(5)

where

t^{'} \in [t_{1}, t_{1} + T]

1 / β

is approximately the number of speckles per pixel (Here

β = 1

Suppose we obtain infinite temporal realizations of

{X (t^{'}), t^{'} \in [t_{1}, t_{1} + T]}

, then we have:

\begin{array}{l} {〈 Y^{2} (t_{1}) 〉}_{e n s e m b l e} = {〈 \frac{1}{T^{2}} \int_{t_{1}}^{t_{1} + T} \int_{t_{1}}^{t_{1} + T} X (t^{'}) X (t^{″}) d t^{'} d t^{″} 〉}_{e n s e m b l e} \\ \approx \frac{1}{T^{2}} {〈 \int_{t_{1}}^{t_{1} + T} \int_{t_{1}}^{t_{1} + T} {〈 X (t^{'}) X (t^{″}) 〉}_{t e m p o r a l} d t^{'} d t^{″} 〉}_{e n s e m b l e} \\ = \frac{1}{T^{2}} \int_{t_{1}}^{t_{1} + T} \int_{t_{1}}^{t_{1} + T} (1 + β {| g_{1} (t^{″} - t^{'}) |}^{2}) d t^{'} d t^{″} {〈 {〈 X (t_{1}) 〉}_{e n s e m b l e}^{2} 〉}_{e n s e m b l e} \\ = \frac{1}{T^{2}} (T^{2} + \int_{0}^{T} \int_{0}^{T} β {| g_{1} (t^{″} - t^{'}) |}^{2} d t^{'} d t^{″}) {〈 μ_{X}^{2} (t_{1}) 〉}_{e n s e m b l e} \\ = \frac{1}{T^{2}} (T^{2} + \int_{0}^{T} \int_{0}^{T} β {| g_{1} (t^{″} - t^{'}) |}^{2} d t^{'} d t^{″}) μ_{Y}^{2} (t_{1}) \\ = {〈 Y (t_{1}) 〉}_{e n s e m b l e}^{2} (1 + \frac{β}{T} \int_{0}^{T} 2 (1 - \frac{τ}{T}) g_{1}^{2} (τ) d τ) \end{array}

(6)

Notes: (i) the temporal average

{〈 X (t^{'}) X (t^{″}) 〉}_{t e m p o r a l}

in the time interval

[t_{1}, t_{1} + T]

is used instead of the function

X (t^{'}) X (t^{″})

in the second step based on the WSS property of

{X (t^{'}), t^{'} \in [t_{1}, t_{1} + T]}

; (ii) Eq. (5) is applied in the third step; (iii) Eq. (3) is used in the fifth step; (iv)

\frac{1}{T^{2}} \int_{0}^{T} \int_{0}^{T} β {| g_{1} (t^{″} - t^{'}) |}^{2} d t^{'} d t^{″} = \frac{β}{T} \int_{0}^{T} 2 (1 - \frac{τ}{T}) g_{1}^{2} (τ) d τ

[15

R. Bandyopadhyay, A. Gittings, S. Suh, P. Dixon, and D. Durian, “Speckle-visibility spectroscopy: A tool to study time-varying dynamics,” Rev. Sci. Instrum. 76(9), 093110 (2005). [CrossRef]

,17

P. K. Dixon and D. J. Durian, “Speckle Visibility Spectroscopy and Variable Granular Fluidization,” Phys. Rev. Lett. 90(18), 184302 (2003). [CrossRef] [PubMed]

] is used in the last step.

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

\frac{{〈 Y^{2} (t_{1}) 〉}_{e n s e m b l e} - {〈 Y (t_{1}) 〉}_{e n s e m b l e}^{2}}{{〈 Y (t_{1}) 〉}_{e n s e m b l e}^{2}} = \frac{β}{T} \int_{0}^{T} 2 (1 - \frac{τ}{T}) g_{1}^{2} (τ) d τ .

(7)

Using

l i m_{N \to + \infty} {〈 Y (t_{1}) 〉}_{e n s e m b l e} \to μ_{Y} (t_{1})

and

l i m_{N \to + \infty} {〈 Y^{2} (t_{1}) 〉}_{e n s e m b l e} - {〈 Y (t_{1}) 〉}_{e n s e m b l e}^{2} \to σ_{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:

K_{Y}^{2} (t_{1}) = \frac{σ_{Y}^{2} (t_{1})}{μ_{Y}^{2} (t_{1})} = \frac{β}{T} \int_{0}^{T} 2 (1 - \frac{τ}{T}) g_{1}^{2} (τ) d τ,

(8)

where

g_{1} (τ)

is the

E (t_{1})

’s autocorrelation function [18

J. W. Goodman, Statistical Optics (Wiley 1985).

Based on different models describing the relation between

v (t)

and

g_{1} (τ)

, the relation between

K_{Y}^{2} (t)

and

v (t)

can also be derived directly [13

A. Fercher and J. Briers, “Flow visualization by means of single-exposure speckle photography,” Opt. Commun. 37(5), 326–330 (1981). [CrossRef]

,15

R. Bandyopadhyay, A. Gittings, S. Suh, P. Dixon, and D. Durian, “Speckle-visibility spectroscopy: A tool to study time-varying dynamics,” Rev. Sci. Instrum. 76(9), 093110 (2005). [CrossRef]

,19

P. Zakharov, A. Völker, A. Buck, B. Weber, and F. Scheffold, “Quantitative modeling of laser speckle imaging,” Opt. Lett. 31(23), 3465–3467 (2006). [CrossRef] [PubMed]

,20

A. B. Parthasarathy, W. J. Tom, A. Gopal, X. Zhang, and A. K. Dunn, “Robust flow measurement with multi-exposure speckle imaging,” Opt. Express 16(3), 1975–1989 (2008). [CrossRef] [PubMed]

]. 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

μ_{Y} (t_{1})

and

σ_{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}) 〉}_{s p a t i a l} = \frac{1}{M} \sum_{j = 0}^{M - 1} {\tilde{I}}_{j} t_{1} \approx {〈 Y (t_{1}) 〉}_{e n s e m b l e}

and

{〈 Y^{2} (t_{1}) 〉}_{s p a t i a l} = \frac{1}{M} \sum_{j = 0}^{M - 1} {\tilde{I}}_{j}^{2} (t_{1}) \approx {〈 Y^{2} (t_{1}) 〉}_{e n s e m b l e}

. The estimations of

μ_{Y} (t_{1})

σ_{Y} (t_{1})

and

K_{Y}^{2} (t_{1})

can be obtained as follows:

K_{s p a t i a l}^{2} (t_{1}) = \frac{σ_{s}^{2} (t_{1})}{μ_{s}^{2} (t_{1})} \approx K_{Y}^{2} (t_{1})

(9)

where

μ_{s} (t_{1}) = {〈 Y (t_{1}) 〉}_{s p a t i a l} \approx μ_{Y} (t_{1})

σ_{s}^{2} (t_{1}) = {〈 Y^{2} (t_{1}) 〉}_{s p a t i a l} - {〈 Y (t_{1}) 〉}_{s p a t i a l}^{2} \approx σ_{Y}^{2} (t_{1})

It is noted that this Eq. (9) is actually the spatial LASCA method [7

J. Briers and S. Webster, “Laser speckle contrast analysis (LASCA): a nonscanning, full-field technique for monitoring capillary blood flow,” J. Biomed. Opt. 1(2), 174–179 (1996). [CrossRef]

]. 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

μ_{s} (t_{1})

and

σ_{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_{μ_{s}} (t_{1})

and

N A_{σ_{s}} (t_{1})

; and (ii) the spatial inhomogeneity of velocities within the window, which could result in estimation noises

N B_{μ_{s}} (t_{1})

and

N B_{σ_{s}} (t_{1})

. Therefore,

μ_{s} (t_{1})

and

σ_{s} (t_{1})

can be modeled as follows:

μ_{s} (t_{1}) = μ_{Y} (t_{1}) + N A_{μ_{s}} (t_{1}) + N B_{μ_{s}} (t_{1}),

(10)

σ_{s} (t_{1}) = σ_{Y} (t_{1}) + N A_{σ_{s}} (t_{1}) + N B_{σ_{s}} (t_{1}) .

(11)

Although a larger window (e.g.

5 \times 5

or larger) can result in smaller

N A_{μ_{s}}

and

N A_{σ_{s}}

based on central limit theory, the

N B_{μ_{s}}

and

N B_{σ_{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_{μ_{s}}

and

N B_{σ_{s}}

, but larger

N A_{μ_{s}}

and

N A_{σ_{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^{'}), t^{'} \in [t_{1}, t_{2}]}

are the same. In practice, we use a CCD camera to record the temporal samples (

{\tilde{I} (t^{'})}

) of integrated intensity random process

Y (t^{'})

\tilde{I} (t_{1} + n Δ t) = \frac{1}{T} \int_{t_{1} + n Δ t}^{t_{1} + n Δ t + T} I (t^{'}) d t^{'},

(12)

where

n = 0, \dots, N - 1

Δ t

is the time interval between two frames. Temporal samples

{\tilde{I} (t^{'} + n Δ t), n = 0, \dots, N - 1}

are thus considered as ensemble samples of

Y (t_{1}) : {〈 Y (t_{1}) 〉}_{t e m p o r a l} = \frac{1}{N} \sum_{n = 0}^{N - 1} \tilde{I} (t_{1} + n Δ t) \approx {〈 Y (t_{1}) 〉}_{e n s e m b l e}

and

{〈 Y^{2} (t_{1}) 〉}_{t e m p o r a l} = \frac{1}{M} \sum_{n = 0}^{M - 1} {\tilde{I}}^{2} (t_{1} + n Δ t) \approx {〈 Y^{2} (t_{1}) 〉}_{e n s e m b l e}

. Then the estimation of

K_{Y}^{2} (t_{1})

, i.e.

μ_{Y} (t_{1})

and

σ_{Y} (t_{1})

individually, is obtained by the temporal estimator defined as follows:

K_{t e m p o r a l}^{2} (t_{1}) = \frac{σ_{t}^{2} (t_{1})}{μ_{t}^{2} (t_{1})} \approx K_{Y}^{2} (t_{1}),

(13)

where

μ_{t} (t_{1}) = {〈 Y (t_{1}) 〉}_{t e m p o r a l} \approx μ_{Y} (t_{1})

σ_{t}^{2} (t_{1}) = {〈 Y^{2} (t_{1}) 〉}_{t e m p o r a l} - {〈 Y (t_{1}) 〉}_{t e m p o r a l}^{2} \approx σ_{Y}^{2} (t_{1})

Again the temporal estimator

K_{t e m p o r a l}^{2}

is actually the temporal LASCA method [8

H. Cheng, Q. Luo, S. Zeng, S. Chen, J. Cen, and H. Gong, “Modified laser speckle imaging method with improved spatial resolution,” J. Biomed. Opt. 8(3), 559–564 (2003). [CrossRef] [PubMed]

]. The temporal LASCA is also contaminated by the estimation noise, i.e. the noise components

N A

and

N B

μ_{t} (t_{1})

and

σ_{t} (t_{1})

as follows:

μ_{t} (t_{1}) = μ_{Y} (t_{1}) + N A_{μ_{t}} (t_{1}) + N B_{μ_{t}} (t_{1}),

(14)

σ_{t} (t_{1}) = σ_{Y} (t_{1}) + N A_{σ_{t}} (t_{1}) + N B_{σ_{t}} (t_{1}) .

(15)

where

N A_{μ t}

and

N A_{σ t}

correspond to the noise due to limited statistics,

N B_{μ t}

and

N B_{σ t}

correspond to the noise due to temporal inhomogeneity of velocities.

Temporal LASCA preserves the spatial resolution because

μ_{t} (t_{1})

and

σ_{t} (t_{1})

are calculated based on

\tilde{I} (t_{1} + n Δ 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

N B

can be ignored. However, if the blood velocity is changing, e.g. elicited in response to functional stimulation, the

N B

noise will be prominent.

To achieve a high performance of statistical calculation (low

N A

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

10 f p s

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

N A

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.

~ 5 m s

). For an interesting pixel, we calculate the

{σ_{s} (t_{1} + n Δ 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_{σ_{s}}

and small

N B_{σ_{s}}

{σ_{s} (t_{1} + n Δ t), n = 0, \dots, N - 1}

. Then, Eq. (11) is simplified into:

σ_{s} (t_{1} + n Δ t) \approx σ_{Y} (t_{1} + n Δ t) + N A_{σ_{s}} (t_{1} + n Δ t)

, where

n = 0, \dots, N - 1

It is noted that each

σ_{Y} (t_{1} + n Δ t)

is related to the velocity

v (t_{1} + n Δ t)

and thus could change during the imaging procedure. The noise component

N A_{σ_{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_{σ_{s}}

can be hypothesized to follow a zero-mean Gaussian distribution, i.e.

N (0, σ_{1}^{2})

, which will be confirmed based on the simulation data in Discussion. Therefore, the discrete sequences

{N A_{σ_{s}} (t_{1} + n Δ 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

{σ_{s} (t_{1} + n Δ t)}

as the random process estimator for averaged

{\bar{σ}}_{Y} (t_{1})

σ_{r p e} (t_{1}) = \frac{1}{N} \sum_{n = 0}^{N - 1} σ_{s} (t_{1} + n Δ t) = {\bar{σ}}_{Y} (t_{1}) + \frac{1}{N} \sum_{n = 0}^{N - 1} N A_{σ_{s}} (t_{1} + n Δ t),

(16)

where

{\bar{σ}}_{Y} (t_{1}) = \frac{1}{N} \sum_{n = 0}^{N - 1} σ_{Y} (t_{1} + n Δ t)

corresponds to the averaged velocity

\bar{v}

during

[t_{1}, t_{1} + (N - 1) Δ t]

. Because

{N A_{σ_{s}} (t_{1} + n Δ t)}

are IID random variables, based on the law of large number,

\lim_{N \to + \infty} \frac{1}{N} \sum_{n = 0}^{N - 1} N A_{σ_{s}} (t_{1} + n Δ t) \to 0

. So, there is a de-noising effect when N is greater than one.

μ_{t} (t_{1})

is used to estimate the averaged

{\bar{μ}}_{Y} (t_{1})

in the random process estimator, so there is no loss of spatial resolution in the calculation of

μ_{r p e} (t_{1})

. Actually, each

\tilde{I} (t_{1} + n Δ t)

can be considered as a ‘spatial estimation’ of

μ_{s} (t_{1} + n Δ t)

with a

1 \times 1

spatial window, which leads to a large noise component

N A_{μ_{s}}

without

N B_{μ_{s}}

noise. Ignoring the

N B

noise in Eq. (10), we have:

\tilde{I} (t_{1} + n Δ t) = μ_{Y} (t_{1} + n Δ t) + N A_{μ_{s}} (t_{1} + n Δ t)

. Then we get:

μ_{r p e} (t_{1}) = \frac{1}{N} \sum_{n = 0}^{N - 1} \tilde{I} (t_{1} + n Δ t) = \frac{1}{N} \sum_{n = 0}^{N - 1} μ_{Y} (t_{1} + n Δ t) + \frac{1}{N} \sum_{n = 0}^{N - 1} N A_{μ_{s}} (t_{1} + n Δ t) .

(17)

Similarly, we model

N A_{μ_{s}} (t_{1} + n Δ t)

as a sequence of IID random variables with a zero-mean Gaussian distribution based on the central limit theorem, i.e.

N (0, σ_{2}^{2})

(this assumption will also be confirmed based on the simulation data in Discussion). Therefore, when

N \to + \infty

\frac{1}{N} \sum_{n = 0}^{N - 1} N A_{μ_{s}} (t_{1} + n Δ t) \to 0

and then

μ_{r p e} (t_{1}) = {\bar{μ}}_{Y} (t_{1}) = \frac{1}{N} \sum_{n = 0}^{N - 1} μ_{Y} (t_{1} + n Δ t)

, where

{\bar{μ}}_{Y} (t_{1})

corresponds to the averaged velocity

\bar{v}

within

[t_{1}, t_{1} + (N - 1) Δ t]

Finally, the new random process estimator is used to estimate

K_{Y}^{2} (t_{1})

for the averaged velocity

\bar{v}

during time interval

[t_{1}, t_{1} + (N - 1) Δ t]

as follows:

K_{r p e}^{2} (t_{1}) = \frac{σ_{r p e}^{2} (t_{1})}{μ_{r p e}^{2} (t_{1})} \approx \frac{{\bar{σ}}_{Y}^{2} (t_{1})}{{\bar{μ}}_{Y}^{2} (t_{1})} = K_{Y}^{2} (t_{1}) .

(18)

Compared with spatial LASCA,

K_{r p e}^{2}

provides much higher spatial resolution because the calculation of

σ_{r p e}

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

N A

and

N B

. 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

N A

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

N B

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

632 n m

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

60 m m

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

5 m s

. Frame rate was

10 f p s

. The resolution of image recorded by the CCD camera was

704 \times 704

pixels corresponding to an imaging area about

4.7 \times 4.7 m 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

2 m m / s

. When the flow was steady after

~ 5 m i n

, 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 (

~ 320 g

) 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

90 m g / k g

of Ketamine and

10 m g / k g

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

70 %

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

5 m m \times 5 m m

area centered at

2.5 m m

lateral and

2.5 m m

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^{°} 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

2 L / m i n

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.5 m A

magnitude and

200 μ s

duration at a frequency of

3 H z

were used for hind paw stimulation. Each stimulation trial consisted of pre-stimulus 200 frames (

20 s

), 200 stimulation frames (

20 s

), and 400 post-stimulus frames (

40 s

), 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

2 m m / s

using the syringe pump. According to Eq. (16) and Eq. (17), the noise parts

\frac{1}{N} \sum_{n = 0}^{N - 1} N A_{μ_{s}} (t_{1} + n Δ t)

and

\frac{1}{N} \sum_{n = 0}^{N - 1} N A_{σ_{s}} (t_{1} + n Δ t)

in the random process estimator would be very closed to zeros with

N = 600

. Therefore, for any pixel,

μ_{r p e}

and

σ_{r p e}

based on all 600 raw images were considered as the true

μ_{Y}

and

σ_{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

σ_{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.

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, 59 s

) while the results of temporal LASCA and random process estimator are calculated using a

1 s

temporal window (10 frames) starting from each second. Figure 1(b) shows the estimations of

σ_{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

N A

(small temporal window) in Eq. (15). Although both LASCA methods are contaminated by

N A

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

N A

noise, there is also

N B

noise in the result of both LASCA methods. Compared to the conventional LASCA methods, random process estimator provides the best estimation of

σ_{Y}

in the sense of smallest fluctuations. Since the estimations of

μ_{Y}

are the same in temporal LASCA and random process estimator, we did not show the results of

μ_{Y}

Figure 1(d) shows the estimations of

K_{Y}^{2}

at each second (

t_{1} = 0, 1, \dots, 59 s

) 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:

S N R = {(\frac{A_{s i g n a l}}{A_{n o i s e}})}^{2} = {(\frac{K_{Y}^{2}}{{[\frac{1}{T} \sum_{n = 0}^{T - 1} {(K_{e s t i m a t i o n}^{2} (n) - K_{Y}^{2})}^{2}]}^{1 / 2}})}^{2},

(19)

where T=60 in this study. For the selected pixel,

K_{Y}^{2} = 0.0053

, the SNRs are

2.75

(temporal LASCA),

10.34

(spatial LASCA), and

31.09

(random process estimator) respectively. For all

60 \times 100

pixels in the white box area in Fig. 1(c), the averaged SNR of random process estimator is

33.64 \pm 6.87

(

m e a n \pm s . t . d .

) while the spatial LASCA and temporal LASCA are only with the averaged SNR

9.08 \pm 2.85

and

3.83 \pm 1.05

respectively. Clearly, random process estimator provides higher denoising performance than either spatial LASCA or temporal LASCA.

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 (

60 s

) 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.

According to Eq. (15), when more frames are used, the noise component

N A

decreases but still exists in the result of temporal LASCA. Furthermore, if there are small changes of blood velocity during the imaging procedure, the

N B

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.

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 (

80 s

) continuously, i.e. 200 pre-stimulus frames (

20 s

), 200 stimulation frames (

20 s

) and 400 post-stimulus frames (

40 s

). 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

F_{n} = \frac{K_{- 20}^{2} - K_{n}^{2}}{K_{- 20}^{2}} \times 100 %,

(20)

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.

E r r o r_{F_{n}} = {(\frac{1}{N_{r o w} N_{c o l}} \sum_{x = 1}^{N_{r o w}} \sum_{y = 1}^{N_{c o l}} F_{n}^{2} (x, y))}^{1 / 2},

(21)

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

(

m e a n \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.

{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 (

{{\bar{F}}_{0}, \dots, {\bar{F}}_{59}}

) is obtained. Figure 6 and Fig. 7 show the averaged changing map

{{\bar{F}}_{0}, \dots, {\bar{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.

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

In the functional imaging experiment, the blood flow in the somatosensory cortex changes in response to electrical stimulation. Therefore, besides the

N A

noise, the

N B

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

4 s

and reaches maximum at

20 s ~ 22 s

, then recovers to the baseline at

32 s

. 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, 59 s}

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, 59 s}

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.

{\bar{F}}_{n}, n = - 19, \dots, 59 s

, 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.

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

σ_{Y}

can be obtained using a

3 \times 3

spatial window. The noise component

N A_{σ_{s}}

of each estimation is then calculated by subtracting the true value

σ_{Y}

. The probability density function (pdf) is then calculated from all 600

N A_{σ_{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}

(

m e a n \pm s . t . d .

). The P-P plot of all 600

N A_{σ_{s}}

in Fig. 9(b) is also very close to that of a Gaussian distribution. The hypothesis that the distribution of

N A_{σ_{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_{σ_{s}}

; (b) the P-P plot of all 600

N A_{σ_{s}}

s; (c) the pdf of

N A_{μ_{s}}

; (d) the P-P plot of all 600

N A_{μ_{s}}

Similarly, the distribution of the noise component

N A_{μ_{s}}

, i.e.

{\tilde{I}}_{i} - μ_{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_{μ_{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

H. Cheng, Y. Yan, and T. Q. Duong, “Laser speckle imaging of rat retinal blood flow with hybrid temporal and spatial analysis method,” Proc. SPIE 7163, 716304 (2009). [CrossRef]

]. Similar to spatial and temporal LASCA, the hybrid method works well when blood flow is constant within both spatial and temporal window, i.e.

N B

noise does not exist. We further process the same data set in our study using the algorithm in [21

H. Cheng, Y. Yan, and T. Q. Duong, “Laser speckle imaging of rat retinal blood flow with hybrid temporal and spatial analysis method,” Proc. SPIE 7163, 716304 (2009). [CrossRef]

] 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

N B

noise while suppressing the

N A

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

N B

noise when there are significant blood flow changes within the temporal window. For example,

N B

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

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.

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.

Acknowledgments

This work is supported by NIH/NIA1R01AG029681. S. Tong is also supported by the New Century Talent Program by the Ministry of Education of China, and Shanghai Shuguang Program (07SG13). P. Miao is also supported by the China Scholarship Council.

References and links

1.	J. Briers, “Laser Doppler, speckle and related techniques for blood perfusion mapping and imaging,” Physiol. Meas. 22(4), 35–66 (2001). [CrossRef]
2.	A. K. Dunn, H. Bolay, M. A. Moskowitz, and D. A. Boas, “Dynamic Imaging of Cerebral Blood Flow Using Laser Speckle,” J. Cereb. Blood Flow Metab. 21(3), 195–201 (2001). [CrossRef] [PubMed]
3.	D. Zhu, W. Lu, Y. Weng, H. Cui, and Q. Luo, “Monitoring thermal-induced changes in tumor blood flow and microvessels with laser speckle contrast imaging,” Appl. Opt. 46(10), 1911–1917 (2007). [CrossRef] [PubMed]
4.	H. Cheng, Y. Yan, and T. Duong, “Temporal statistical analysis of laser speckle images and its application to retinal blood-flow imaging,” Opt. Express 16(14), 214–219 (2008). [CrossRef]
5.	A. Kharlamov, B. R. Brown, K. A. Easley, and S. C. Jones, “Heterogeneous response of cerebral blood flow to hypotension demonstrated by laser speckle imaging flowmetry in rats,” Neurosci. Lett. 368(2), 151–156 (2004). [CrossRef] [PubMed]
6.	T. Durduran, M. G. Burnett, G. Yu, C. Zhou, D. Furuya, A. G. Yodh, J. A. Detre, and J. H. Greenberg, “Spatiotemporal Quantification of Cerebral Blood Flow During Functional Activation in Rat Somatosensory Cortex Using Laser-Speckle Flowmetry,” J. Cereb. Blood Flow Metab. 24(5), 518–525 (2004). [CrossRef] [PubMed]
7.	J. Briers and S. Webster, “Laser speckle contrast analysis (LASCA): a nonscanning, full-field technique for monitoring capillary blood flow,” J. Biomed. Opt. 1(2), 174–179 (1996). [CrossRef]
8.	H. Cheng, Q. Luo, S. Zeng, S. Chen, J. Cen, and H. Gong, “Modified laser speckle imaging method with improved spatial resolution,” J. Biomed. Opt. 8(3), 559–564 (2003). [CrossRef] [PubMed]
9.	P. Miao, M. Li, H. Fontenelle, A. Bezerianos, Y. Qiu, and S. Tong, “Imaging the Cerebral Blood Flow with Enhanced Laser Speckle Contrast Analysis (eLASCA) by Monotonic Point Transformation,” IEEE Trans. Biomed. Eng. 56(4), 1127–1133 (2009). [CrossRef] [PubMed]
10.	P. Lemieux and D. Durian, “Investigating non-Gaussian scattering processes by using nth-order intensity correlation functions,” J. Opt. Soc. Am. A 16(7), 1651–1664 (1999). [CrossRef]
11.	D. Boas and A. Yodh, “Spatially varying dynamical properties of turbid media probed with diffusing temporal light correlation,” J. Opt. Soc. Am. A 14(1), 192–215 (1997). [CrossRef]
12.	B. Berne, and R. Pecora, Dynamic Light Scattering: with Applications to Chemistry, Biology and Physics (Dover Publications, 2000).
13.	A. Fercher and J. Briers, “Flow visualization by means of single-exposure speckle photography,” Opt. Commun. 37(5), 326–330 (1981). [CrossRef]
14.	D. J. Pine, D. A. Weitz, P. M. Chaikin, and E. Herbolzheimer, “Diffusing wave spectroscopy,” Phys. Rev. Lett. 60(12), 1134–1137 (1988). [CrossRef] [PubMed]
15.	R. Bandyopadhyay, A. Gittings, S. Suh, P. Dixon, and D. Durian, “Speckle-visibility spectroscopy: A tool to study time-varying dynamics,” Rev. Sci. Instrum. 76(9), 093110 (2005). [CrossRef]
16.	G. Grimmett, and D. Stirzaker, Probability and random processes (Oxford University Press, 2001).
17.	P. K. Dixon and D. J. Durian, “Speckle Visibility Spectroscopy and Variable Granular Fluidization,” Phys. Rev. Lett. 90(18), 184302 (2003). [CrossRef] [PubMed]
18.	J. W. Goodman, Statistical Optics (Wiley 1985).
19.	P. Zakharov, A. Völker, A. Buck, B. Weber, and F. Scheffold, “Quantitative modeling of laser speckle imaging,” Opt. Lett. 31(23), 3465–3467 (2006). [CrossRef] [PubMed]
20.	A. B. Parthasarathy, W. J. Tom, A. Gopal, X. Zhang, and A. K. Dunn, “Robust flow measurement with multi-exposure speckle imaging,” Opt. Express 16(3), 1975–1989 (2008). [CrossRef] [PubMed]
21.	H. Cheng, Y. Yan, and T. Q. Duong, “Laser speckle imaging of rat retinal blood flow with hybrid temporal and spatial analysis method,” Proc. SPIE 7163, 716304 (2009). [CrossRef]

OCIS Codes
(110.4280) Imaging systems : Noise in imaging systems
(110.6150) Imaging systems : Speckle imaging
(170.3880) Medical optics and biotechnology : Medical and biological imaging

ToC Category:
Medical Optics and Biotechnology

History
Original Manuscript: October 26, 2009
Revised Manuscript: December 7, 2009
Manuscript Accepted: December 14, 2009
Published: December 23, 2009

Virtual Issues
Vol. 5, Iss. 2 Virtual Journal for Biomedical Optics

Citation
Peng Miao, Nan Li, Nitish V. Thakor, and Shanbao Tong, "Random process estimator for laser speckle imaging of cerebral blood flow," Opt. Express 18, 218-236 (2010)
http://www.opticsinfobase.org/vjbo/abstract.cfm?URI=oe-18-1-218

Sort: Author | Year | Journal | Reset

References

J. Briers, “Laser Doppler, speckle and related techniques for blood perfusion mapping and imaging,” Physiol. Meas. 22(4), 35–66 (2001). [CrossRef]
A. K. Dunn, H. Bolay, M. A. Moskowitz, and D. A. Boas, “Dynamic Imaging of Cerebral Blood Flow Using Laser Speckle,” J. Cereb. Blood Flow Metab. 21(3), 195–201 (2001). [CrossRef] [PubMed]
D. Zhu, W. Lu, Y. Weng, H. Cui, and Q. Luo, “Monitoring thermal-induced changes in tumor blood flow and microvessels with laser speckle contrast imaging,” Appl. Opt. 46(10), 1911–1917 (2007). [CrossRef] [PubMed]
H. Cheng, Y. Yan, and T. Duong, “Temporal statistical analysis of laser speckle images and its application to retinal blood-flow imaging,” Opt. Express 16(14), 214–219 (2008). [CrossRef]
A. Kharlamov, B. R. Brown, K. A. Easley, and S. C. Jones, “Heterogeneous response of cerebral blood flow to hypotension demonstrated by laser speckle imaging flowmetry in rats,” Neurosci. Lett. 368(2), 151–156 (2004). [CrossRef] [PubMed]
T. Durduran, M. G. Burnett, G. Yu, C. Zhou, D. Furuya, A. G. Yodh, J. A. Detre, and J. H. Greenberg, “Spatiotemporal Quantification of Cerebral Blood Flow During Functional Activation in Rat Somatosensory Cortex Using Laser-Speckle Flowmetry,” J. Cereb. Blood Flow Metab. 24(5), 518–525 (2004). [CrossRef] [PubMed]
J. Briers and S. Webster, “Laser speckle contrast analysis (LASCA): a nonscanning, full-field technique for monitoring capillary blood flow,” J. Biomed. Opt. 1(2), 174–179 (1996). [CrossRef]
H. Cheng, Q. Luo, S. Zeng, S. Chen, J. Cen, and H. Gong, “Modified laser speckle imaging method with improved spatial resolution,” J. Biomed. Opt. 8(3), 559–564 (2003). [CrossRef] [PubMed]
P. Miao, M. Li, H. Fontenelle, A. Bezerianos, Y. Qiu, and S. Tong, “Imaging the Cerebral Blood Flow with Enhanced Laser Speckle Contrast Analysis (eLASCA) by Monotonic Point Transformation,” IEEE Trans. Biomed. Eng. 56(4), 1127–1133 (2009). [CrossRef] [PubMed]
P. Lemieux and D. Durian, “Investigating non-Gaussian scattering processes by using nth-order intensity correlation functions,” J. Opt. Soc. Am. A 16(7), 1651–1664 (1999). [CrossRef]
D. Boas and A. Yodh, “Spatially varying dynamical properties of turbid media probed with diffusing temporal light correlation,” J. Opt. Soc. Am. A 14(1), 192–215 (1997). [CrossRef]
B. Berne, and R. Pecora, Dynamic Light Scattering: with Applications to Chemistry, Biology and Physics (Dover Publications, 2000).
A. Fercher and J. Briers, “Flow visualization by means of single-exposure speckle photography,” Opt. Commun. 37(5), 326–330 (1981). [CrossRef]
D. J. Pine, D. A. Weitz, P. M. Chaikin, and E. Herbolzheimer, “Diffusing wave spectroscopy,” Phys. Rev. Lett. 60(12), 1134–1137 (1988). [CrossRef] [PubMed]
R. Bandyopadhyay, A. Gittings, S. Suh, P. Dixon, and D. Durian, “Speckle-visibility spectroscopy: A tool to study time-varying dynamics,” Rev. Sci. Instrum. 76(9), 093110 (2005). [CrossRef]
G. Grimmett, and D. Stirzaker, Probability and random processes (Oxford University Press, 2001).
P. K. Dixon and D. J. Durian, “Speckle Visibility Spectroscopy and Variable Granular Fluidization,” Phys. Rev. Lett. 90(18), 184302 (2003). [CrossRef] [PubMed]
J. W. Goodman, Statistical Optics (Wiley 1985).
P. Zakharov, A. Völker, A. Buck, B. Weber, and F. Scheffold, “Quantitative modeling of laser speckle imaging,” Opt. Lett. 31(23), 3465–3467 (2006). [CrossRef] [PubMed]
A. B. Parthasarathy, W. J. Tom, A. Gopal, X. Zhang, and A. K. Dunn, “Robust flow measurement with multi-exposure speckle imaging,” Opt. Express 16(3), 1975–1989 (2008). [CrossRef] [PubMed]
H. Cheng, Y. Yan, and T. Q. Duong, “Laser speckle imaging of rat retinal blood flow with hybrid temporal and spatial analysis method,” Proc. SPIE 7163, 716304 (2009). [CrossRef]

Cited By

Alert me when this paper is cited

OSA is able to provide readers links to articles that cite this paper by participating in CrossRef's Cited-By Linking service. CrossRef includes content from more than 3000 publishers and societies. In addition to listing OSA journal articles that cite this paper, citing articles from other participating publishers will also be listed.

Click here to see a list of articles that cite this paper