Robust quantitative single-exposure laser speckle imaging with true flow speckle contrast in the temporal and spatial domains

Chenge Wang; Zili Cao; Xin Jin; Weihao Lin; Yang Zheng; Bixin Zeng; M. Xu

doi:10.1364/BOE.10.004097

1. Introduction

When coherent light interacts with a turbid medium, the interference between the outgoing waves produces grainy speckle patterns which encode the phase fluctuation of all rays (random phasors) reaching a point. The contrast of laser speckles reduces with the motion of scatterers inside the turbid medium. Laser speckle contrast hence can be used to infer the dynamic property of the medium. Laser speckle contrast imaging (LSCI, see recent reviews [1–4]) has now been widely used in monitoring blood flow in brain, skin, retina, arthrosis and etc due to advantages including simplicity, high spatial and temporal resolution, and large field of view without scanning [5–8].

Although LSCI has a wide range of applications and a long history, the recovery of absolute flow velocity from LSCI measurements remains a challenge, especially when the measurement is compounded by static scattering and noise. For static scattering in laser speckle imaging, Li et al. showed that the static scattering effect can be partially suppressed by using the temporal rather than spatial contrast analysis of laser speckles [9] as the static scattering is an invariant quantity with time. Zakharov et al. [10,11] presented a data processing scheme to correctly separate dynamic and static components within the speckle contrast based on their different decorrelation behavior across speckle patterns captured at consecutive times. Parthasarathy et al. [12,13] and Kazmi et al. [6,14] later demonstrated a multi-exposure laser speckle contrast imaging method, which quantifies and eliminates the influence of static scattering from speckle contrasts measured under different exposure times using a laser speckle contrast model. For LSCI measurement noise, the correction of the variance of the shot noise and sensor dark currents were found to be crucial to estimate the true speckle contrast [15,16]. Yuan et al. [16] increased the signal-to-noise ratio (SNR) of LSCI with noise correction to detect small blood flow changes caused by brain activity. A systematic study and recommended practical recipe to obtain true flow velocity from LSCI measurements addressing both static scattering and measurement noise is, however, still lacking.

In this article, we analysed laser speckle flow imaging from the first principle and provided a complete procedure covering the LSCI system calibration, static scattering removal, and measurement noise estimation and correction to obtain a genuine flow speckle contrast and the flow speed from single-exposure LSCI measurements. We demonstrated the power of our LSCI analysis recipe by imaging phantom flow at varying speeds. Experimental results show that our procedure greatly enhances the linear sensitivity of the flow index (defined as the inverse decorrelation time) and the linearity covers the full span of flow speeds from 0 mm/s to 40 mm/s. The true flow speed is recovered regardless of the overlying static scattering layers and the type of statistics (temporal or spatial). The proposed LSCI analysis framework hence paves the way to estimate the true flow speed in the wide array of laser speckle contrast imaging applications.

2. Theory and Data Analysis

2.1 Theoretical basis

The spatial intensity distribution of the speckle pattern fluctuates with the motion of the scattering particles under the illumination of coherent light. The recorded pattern by a camera is the integration of all instantaneous speckles over the exposure time. The faster the scattering particles move, the more blurred the recorded pattern becomes. The degree of blurring is quantified by the contrast [17] given by

K = \frac{\sqrt{\bar{\sum {(I - \bar{I})}^{2}}}}{\bar{I}}

where

\bar{I}

represents the mean of light intensity I over a small region (spatial contrast) or over a short durance of time (temporal contrast). For “fully developed” static speckles, the spatial contrast K equals to 1.

We will assume the scattered electric field containing both dynamic and static components

E (t) e^{- i ω t} = [E_{f} (t) + E_{s}] e^{- i ω t}

with ω being the angular frequency of light. The dynamic component consists of photons which have at least been scattered by moving scatterers (flow) once and the static component consists of photons being scattered by static scatterers alone.

The electric field temporal autocorrelation function can be written as

〈 E (t) E^{*} (t + τ) 〉 = G_{1} (τ) + I_{S}

where

〈 〉

means average over t,

G_{1} (τ) = 〈 E_{f} (t) E_{f}^{*} (t + τ) 〉

is the electric field temporal autocorrelation function related to flow, and

I_{s} = {| E_{s} |}^{2}

.

In practice only light intensity fluctuation signals can be recorded. The intensity autocorrelation function is defined as $G_{2} (τ) = 〈 I (t) I (t + τ) 〉$ where $I (t) = I_{f} (t) + I_{s}$ and $I_{f} (t) = {| E_{f} (t) |}^{2}$ assuming the dynamic and static electric fields are uncorrelated, i.e., $〈 E_{f} (t) E_{s}^{*} 〉 = 0$ . In terms of the normalized dynamic electric field and full intensity autocorrelations $g_{1} (τ) = G_{1} (τ) / \bar{I_{f}}$ and $g_{2} (τ) = G_{2} (τ) / G_{2} (0)$ , the full Siegert relation [18,19] expresses

g_{2} (τ) = 1 + \frac{β}{{({\bar{I}}_{f} + I_{s})}^{2}} [{\bar{I}}_{f}^{2} {| g_{1} (τ) |}^{2} + 2 {\bar{I}}_{f} I_{s} | g_{1} (τ) |]

where

{\bar{I}}_{f} = 〈 I_{f} (t) 〉

,

β \leq 1

is a parameter that accounts for the reduction in the measured contrast due to averaging (by the detector) over uncorrelated speckles. Note

g_{1}

is real and non-negative [19].

The speckle contrast under an exposure time T, expressed as $K^{2} = T^{- 2} 〈 \int_{0}^{T} \int_{0}^{T} I (t_{1}) I (t_{2}) d t_{1} d t_{2} 〉 / {〈 I 〉}^{2} - 1$ , reduces now to

K^{2} = ρ^{2} \frac{2 β}{T} \int_{0}^{T} (1 - \frac{τ}{T}) {| g_{1} (τ) |}^{2} d τ + 2 ρ (1 - ρ) \frac{2 β}{T} \int_{0}^{T} (1 - \frac{τ}{T}) | g_{1} (τ) | d τ

where the average intensity

〈 I 〉 = {\bar{I}}_{f} + I_{s}

and the dynamic fraction

ρ = {\bar{I}}_{f} / ({\bar{I}}_{f} + I_{s})

[12]. We will use

K^{2}

as the speckle contrast instead of the conventional contrast K as the former relates simply to the flow velocity and is with additive noise alone.

Equation (5) is the theoretical temporal contrast from the random process taken by the photons migrating through a turbid medium. Some complexities arise when evaluating $K^{2}$ from measurement data. First the measurement noise (of zero mean) introduces an extra variance term $κ_{noise}^{2}$ . Second when using spatial ensemble average for the evaluation of the contrast rather than temporal statistics, the extra terms appearing in $K^{2}$ will be $κ_{noise}^{2} + κ_{ne}^{2}$ with $κ_{ne}^{2}$ being the non-ergodic contribution from the static field. This motivates us to introduce the dynamic (flow) contrast $K_{f}^{2} = 〈 I_{f}^{2} 〉 / {〈 I_{f} 〉}^{2} - 1$ defined in terms of the dynamic component alone. The measured speckle contrast can be expressed as

K^{2} = ρ^{2} K_{f}^{2} + κ_{noise}^{2}

in the temporal domain and

K^{2} = ρ^{2} K_{f}^{2} + κ_{noise}^{2} + κ_{ne}^{2}

in the spatial domain. The value of

κ_{noise}^{2}

and

κ_{ne}^{2}

can be evaluated from calibration and measurement data as shown later. The dynamic contrast remains the same when evaluated in either the temporal or spatial domain.

A velocity distribution model for the moving particles is needed to relate $K_{f}^{2}$ to the flow speed. With the commonly used Lorentz velocity distribution model, the dynamic electric field autocorrelation function can be written as $g_{1} (τ) = \exp (- t / τ_{c})$ [12], yielding

\begin{matrix} K_{f}^{2} = β \frac{\exp (- 2 x) - 1 + 2 x}{2 x^{2}} + 4 β (ρ^{- 1} - 1) \frac{\exp (- x) - 1 + x}{x^{2}} \\ = β (2 ρ^{- 1} - 1) - \frac{2}{3} β ρ^{- 1} x + \frac{1}{6} β (ρ^{- 1} + 1) x^{2}, x ≪ 1 \\ = β (4 ρ^{- 1} - 3) \frac{1}{x} + β (\frac{7}{2} - 4 ρ^{- 1}) \frac{1}{x^{2}}, x ≫ 1 \end{matrix}

where

x = T / τ_{c}

and

τ_{c}

is the decorrelation time. The inverse decorrelation time increases with the flow speed and can be regarded as the flow index.

In the next section, we will examine system calibration, sample measurement and data analysis to provide a complete procedure for static scattering removal, and measurement noise estimation and correction to obtain the true flow speckle contrast $K_{f}^{2}$ and the flow speed from single-exposure LSCI measurements.

2.2 System calibration, measurement, and data analysis

Let's consider a set of speckle images $I_{i}$ ( $1 \leq i \leq N$ ) at time $t_{i} = i Δ t$ with an exposure time T, i.e.,

I_{i} (x, y) = I_{s} (x, y) + I_{f i} (x, y) + n_{i} (x, y)

Here the recorded image consists of the static component

I_{s}

, the dynamic component

I_{f}

and the noise n. The noise [15] presented in the measurement is mainly comprised of the dark counts

n_{d}

and the shot noise

n_{s}

, i.e.,

n_{i} (x, y) = n_{d i} (x, y) + n_{s i} (x, y)

.

The dark counts $n_{d}$ and the variance of dark counts $V a r (n_{d}) = 〈 n_{d}^{2} 〉 - {〈 n_{d} 〉}^{2}$ can be easily acquired by taking multiple dark frames at the same exposure time and camera gain (with all light off). One could use the temporal average to get the dark count and its variance pixel by pixel when the number of the dark frames is large ( $\geq 50$ ) or use the spatial average over a sliding $N_{P} \times N_{P}$ (typical $N_{P}$ = 7) window otherwise. If the behavior of dark counts is assumed uniform across the whole sensor frame, the mean and the variance of the dark counts are given by further averaging over the whole sensor frame. In many cameras, the recorded intensity has been pre-subtracted by certain base. In this case, $〈 n_{d} 〉$ should be estimated by the median value and the variance $V a r (n_{d}) = 2 〈 n_{d}'^{2} 〉$ where $n_{d}' = n_{d} - 〈 n_{d} 〉$ with the negative values of $n_{d}^{}'$ replaced by zeros.

We would always subtract the dark counts from $I_{i}$ before further analysis. After this pre-processing, the speckle image becomes $I_{i}^{c} (x, y) = I_{i} (x, y) - 〈 n_{d} (x, y) 〉$ and the noise term is replaced by $n_{i}^{c} (x, y) = n_{i} (x, y) - 〈 n_{d} (x, y) 〉$ . The noise satisfies $〈 n_{i}^{c} 〉 = 0$ and $V a r (n_{i}^{c}) = V a r (n_{s i}) + V a r (n_{d})$ . The shot noise obeys a Poisson distribution with the mean $〈 n_{s i} 〉 = 0$ and a variance equal to $〈 I_{i}^{c} (x, y) 〉 / γ$ as the camera converts the photoelectrons to digital counts where $γ$ is the analog to digital conversion factor [15]. The $γ$ factor is typically the same across the sensor frame and thus is obtained by further averaging over the sensor frame. Under such a shot noise model,

V a r (n_{i}^{c}) = α {〈 I_{i}^{c} (x, y) 〉}^{2} + 〈 I_{i}^{c} (x, y) 〉 / γ + V a r (n_{d})

where an extra quadratic term in

〈 I_{i}^{c} (x, y) 〉

is added to account for other noise sources such as the laser fluctuations.

The temporal and spatial averages of the sequence of single-exposure speckle images then satisfy:

{〈 I_{i}^{c} (x, y) 〉}_{i} = I_{s} (x, y) + {〈 I_{f i} (x, y) 〉}_{i}

{〈 I_{i}^{c} (x, y) I_{j}^{c} (x, y) 〉}_{i} = I_{s}^{2} (x, y) + {〈 I_{f i} (x, y) I_{f j} (x, y) 〉}_{i} + 2 I_{s} (x, y) {〈 I_{f i} (x, y) 〉}_{i}

and

{〈 I_{i}^{c 2} (x, y) 〉}_{i} = I_{s}^{2} (x, y) + {〈 I_{f i}^{2} (x, y) 〉}_{i} + {〈 n_{i}^{c 2} (x, y) 〉}_{i} + 2 I_{s} (x, y) {〈 I_{f i} (x, y) 〉}_{i}

with temporal statistics and

{〈 I_{i}^{c} (x, y) 〉}_{x y} = {〈 I_{s} (x, y) 〉}_{x y} + {〈 I_{f i} (x, y) 〉}_{x y}

{〈 I_{i}^{c} (x, y) I_{j}^{c} (x, y) 〉}_{x y} = {〈 I_{s}^{2} (x, y) 〉}_{x y} + {〈 I_{f i} (x, y) I_{f j} (x, y) 〉}_{x y} + 2 {〈 I_{s} (x, y) 〉}_{x y} {〈 I_{f i} (x, y) 〉}_{x y}

and

{〈 I_{i}^{2} (x, y) 〉}_{x y} = {〈 I_{s}^{2} (x, y) 〉}_{x y} + {〈 I_{f i}^{2} (x, y) 〉}_{x y} + 2 {〈 I_{s} (x, y) 〉}_{x y} {〈 I_{f i} (x, y) 〉}_{x y} + {〈 n_{i}^{c 2} (x, y) 〉}_{x y}

with spatial statistics. Here

j \equiv i + Δ

with

Δ \neq 0

,

{〈 〉}_{i}

means averaging over the N temporal instances and

{〈 〉}_{x y}

means the spatial average over a sliding

N_{P} \times N_{P}

window. We have used the fact

{〈 I_{f i} (x, y) 〉}_{x y} = {〈 I_{f j} (x, y) 〉}_{x y}

due to the ergodic nature of the dynamic component. When the time difference satisfies

t_{j} - (t_{i} + T) ≫ τ_{c}

, the complete decorrelation between the dynamic component measured at two different times ensures the important identities

{〈 I_{f i} (x, y) I_{f, i + Δ} (x, y) 〉}_{i} = {〈 I_{f i} (x, y) 〉}_{i} {〈 I_{f, i + Δ} (x, y) 〉}_{i}

and

{〈 I_{f i} (x, y) I_{f j} (x, y) 〉}_{x y} = {〈 I_{f i} (x, y) 〉}_{x y} {〈 I_{f j} (x, y) 〉}_{x y}

as

E_{f i} (x, y)

is a zero-mean Gaussian variable. These identities could serve as the data consistency check.

One important consequence is that for any two speckle images $I_{i}$ and $I_{j}$ taken at times satisfying $t_{j} - (t_{i} + T) ≫ τ_{c}$ , the following holds

\frac{{〈 I_{i}^{c} (x, y) I_{j}^{c} (x, y) 〉}_{x y} - {〈 I_{i}^{c} (x, y) 〉}_{x y} {〈 I_{j}^{c} (x, y) 〉}_{x y}}{{〈 I_{i}^{c} (x, y) 〉}_{x y}^{2}} = β (N_{p}) {(1 - ρ)}^{2}

where the coherence factor of the imaging system is defined as

β (N_{p}) = \frac{{〈 I_{s}^{2} (x, y) 〉}_{x y}}{{〈 I_{s} (x, y) 〉}_{x y}^{2}} - 1

associated with spatial averaging over the sliding window of size

N_{P}

. Equation (17) is also correct for a pure static scattering sample producing fully developed speckles (ρ = 0) as long as

j \neq i

.

2.2.1 System calibration

As stated earlier, the dark counts $n_{d}$ and the variance of dark counts $V a r (n_{d})$ is first acquired by taking multiple dark frames at the identical experimental condition (the same exposure time, camera gain and etc) with all light off. The other system parameters (the coherence factor $β$ and the behavior of $V a r (n_{}^{c})$ vs $〈 I_{}^{c} 〉$ ) can be directly evaluated from a set of fully developed speckle images taken on a pure static scattering sample such as a reflection standard. Indeed, according to Eq. (17), we have

β (N_{p}) = \frac{{〈 I_{i}^{c} (x, y) I_{i + Δ}^{c} (x, y) 〉}_{x y}}{{〈 I_{i}^{c} (x, y) 〉}_{x y} {〈 I_{i + Δ}^{c} (x, y) 〉}_{x y}} - 1

for

Δ \neq 0

in this case. Here the spatial average should use the largest window size (full image if possible) due to the reason discussed in [20] for the temporal speckle contrast. For the spatial speckle contrast

N_{P}

in Eq. (19) should take the same value used for the contrast calculation. The

β

factor is typically the same across the sensor frame and thus obtained by further averaging over the sensor frame.

The noise variance $V a r (n_{}^{c})$ associated with a particular $〈 I_{}^{c} 〉$ can be found as

V a r (n^{c}) = {〈 I_{i}^{c 2} (x, y) 〉}_{i} - {〈 I_{i}^{c} (x, y) 〉}_{i}^{2}

or

V a r (n^{c}) = \sqrt{{〈 I_{i}^{c 2} (x, y) 〉}_{x y} {〈 I_{i + Δ}^{c 2} (x, y) 〉}_{x y}} - {〈 I_{i}^{c} (x, y) I_{i + Δ}^{c} (x, y) 〉}_{x y}

with the temporal or spatial statistics, respectively. Multiple sets of such speckle images under the identical experimental condition and varying incident intensities are measured to cover the full range of

〈 I_{}^{c} 〉

. By fitting to Eq. (10), the system noise behavior is then characterized.

We note the above results on $β$ and $V a r (n_{}^{c})$ should be independent of the $Δ \neq 0$ . In reality a slight dependence on ∆ can be observed due to the inevitable system instability. In this case, the correct values of $β$ and $V a r (n_{}^{c})$ are obtained by extrapolating to $Δ = T / 2 Δ t$ .

2.2.2 Sample measurement and data processing

The sample containing both dynamic and static scatterers are then imaged under the identical experimental condition to yield a new set of dark current removed speckle images $I_{i}^{c} (x, y)$ ( $1 \leq i \leq N$ ).

The dynamic fraction $ρ$ can be determined using Eq. (17) as

{(1 - ρ)}^{2} = \frac{1}{β (N_{p})} [\frac{{〈 I_{i}^{c} (x, y) I_{i + Δ}^{c} (x, y) 〉}_{x y}}{{〈 I_{i}^{c} (x, y) 〉}_{x y}^{2}} - 1]

for

t_{i + Δ} - (t_{i} + T) ≫ τ_{c}

.

The noise level can also be estimated directly from the set of speckle images via

κ_{n o i s e}^{2} = \frac{α {〈 I_{i}^{c} (x, y) 〉}_{i}^{2} + {〈 I_{i}^{c} (x, y) 〉}_{i} / γ + V a r (n_{d})}{{〈 I_{i}^{c} (x, y) 〉}_{i}^{2}}

in the temporal domain or

κ_{n o i s e}^{2} = \frac{α {〈 I_{i}^{c} (x, y) 〉}_{x y}^{2} + {〈 I_{i}^{c} (x, y) 〉}_{x y} / γ + V a r (n_{d})}{{〈 I_{i}^{c} (x, y) 〉}_{x y}^{2}}

in the spatial domain (they are equal by ergodicity).

Using Eqs. (13) and (16), the temporal speckle contrast is then

K^{2} = \frac{{〈 I_{i}^{c 2} (x, y) 〉}_{i}}{{〈 I_{i}^{c} (x, y) 〉}_{i}^{2}} - 1 = ρ^{2} K_{f}^{2} + κ_{n o i s e}^{2}

and the spatial speckle contrast is

K^{2} = \frac{{〈 I_{i}^{c 2} (x, y) 〉}_{x y}}{{〈 I_{i}^{c} (x, y) 〉}_{x y}^{2}} - 1 = ρ^{2} K_{f}^{2} + β (N_{p}) {(1 - ρ)}^{2} + κ_{n o i s e}^{2}

from which we can identify the non-ergodic contribution from the static field to be

κ_{ne}^{2} = β (N_{p}) {(1 - ρ)}^{2}

Finally, the velocity information of the sample can be obtained by solving the flow contrast $K_{f}^{2}$ and fitting to Eq. (8) to obtain the decorrelation time $τ_{c}$ and the flow index.

Figure 1 outlines the complete procedure of system calibration, sample measurement, noise correction, and static scattering removal for robust quantitative single-exposure laser speckle imaging.

Fig. 1 Experimental and data analysis framework for robust quantitative single-exposure laser speckle imaging.

Abstract

1. Introduction

2. Theory and Data Analysis

2.1 Theoretical basis

2.2 System calibration, measurement, and data analysis

2.2.1 System calibration

2.2.2 Sample measurement and data processing

2.2.3 Experimental setup

3. Results

3.1 Results of system calibration

3.2 Results of dynamic sample measurements

3.2.1 Importance of static scattering removal and noise correction

3.2.2 Effects of different static scattering

3.2.3 Temporal speckles vs spatial speckles

3.2.4 2D flow profile

4. Discussions

5. Conclusion

Funding

Disclosures

References

Cited By

Figures (16)

Tables (1)

Equations (28)

Biomedical Optics Express