Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Ensemble averaging laser speckle contrast imaging: statistical model of improvement as function of static scatterers

Open Access Open Access

Abstract

The appearance of the common artifacts of laser speckle contrast imaging (LSCI), namely the granularity in flow rate estimation caused by static scatterers, is a well-known phenomenon. This artifact can be greatly reduced in spatial speckle contrast calculation using interframe decorrelated illumination, forcing true ensemble averaging. We propose a statistical model, which describes the effect of multiple image acquisitions on the contrast map quality when the illumination stable and when the illumination is decorrelated frame by frame. We investigate the improvement as a function of the ratio of dynamic and static scatterers by formulating a statistical distribution based model, using in simulation, flow phantom and in vivo experiments. Our main finding is that the ensemble averaging yields limited improvement in several practical cases due to the highly heterogeneous scatterer structure of living tissues.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

1. Introduction

The first description of the Laser Speckle Contrast Imaging (LSCI) was published as early as 1981 [1] regarding the study of the perfusion of blood flow. However, the field of application is wide as it can be used in the medical field for non-invasive measurement of microcirculation and blood perfusion in the brain surface [2], skin [3] and retina [4], as well as estimating semen motility at artificial insemination centers and during fertility tests [5]. In agriculture, dynamic laser speckle or biospeckle is used to monitor the biological activity of the sample, including but not limited to identifying fungi contamination of seeds and fruits [6,7], quantify plant maturation [8], and to monitor fruit and vegetable quality [9]. The transportation of water and nutrients in the venation of leaves [10] can be observed similarly to blood flow in microcirculation.

Static granularity in the contrast map caused by static scatterers is a well-known phenomenon. To overcome this artifact, as active noise reduction technique, interframe decorrelated illumination was introduced in [11], which, in effect, is capable of performing ensemble averaging. Though, the technique is highly cited and promises significant spatial noise reduction, since its introduction, to our best knowledge no actual analysis has been performed to compare the steady illumination and ensemble averaging techniques in biological samples.

Laser speckles are the result of the superposition of scattered coherent light forming interference patterns. The movement of the scatterers creates phase fluctuations, and as a result, the interference patterns decorrelate during observation time. Inert materials possess a behavior of unvarying patterns, while decorrelation will increase with the movement of the scatterers and with other local activities of the living tissue. Among the many analyzing techniques (e.g. Fujii [12], Generalized Difference [13], Absolute Value Difference [14]). In LSCI, to quantify decorrelation and flow rate, the speckle contrast is calculated on a digital image or series of images of the sample. For the calculation, as the basis of the LSCI technology, the standard deviance of the intensity is normalized by the mean intensity in spatial or temporal sliding windows to form a contrast map. The spatial speckle contrast is defined as:

$$K_s = \frac{\sigma(I_s)}{\langle I_s \rangle},$$
where $\sigma (I_s)$ and $\langle I_s \rangle$ are the standard deviation and the sample mean of a small spatial neighborhood of $n_s \times n_s$ intensity values. The reason for the small neighborhood (typically $n_s$ is set to 3 to 9) is to maintain the spatial resolution of the result. Speckle contrast calculation and postprocessing to estimate the flow rate are heavily affected and have high statistical uncertainty due to this small sample size. Thus, it is common practice to capture a sequence of frames ($n \approx 20-100$), the spatial speckle contrast is calculated separately on each of them and then averaged across the contrast maps gathering information of an $n_s \times n_s \times n$ cuboid. Another fundamental approach is to use temporal or combined spatiotemporal contrast calculations [1517]. The temporal methods rely on determining the first-order statistics of the speckle fluctuations data of a sequence [16]:
$$K_t = \frac{\sigma(I_t)}{\langle I_t \rangle} = \sqrt{\frac{\langle I_t^{2} \rangle - \langle I_t \rangle^{2}}{\langle I_t \rangle^{2}}},$$
where $\langle I_t \rangle$ and $\langle I_t^{2} \rangle$ are the mean and mean square values of the speckle intensity variations during exposure time. The temporal contrast calculation works on a $1 \times 1 \times n$ neighborhood using $n$ temporal length that depends on the exposure time and frame rate [16]. The exact relation between the contrast and flow rate is a complex question and depends on multiple factors [1723].

First, we describe the role of ensemble averaging in LSCI, then we provide the theory behind the technique and the effect on static and dynamic scatterers to describe its expected influence on flowmetry. The proposed statistical model is validated on numerical simulations, and in channel slide experiments. The effect of ensemble averaging is demonstrated on three experiments:

  • • microfluidic - Intralipid in a channel slide phantom
  • • in vivo - blood flow in vessel of mouse brain surface
  • • in vivo - flow visualisation on the underside of a fresh Ficus benjamina leaf

2. Ensemble averaging

In case of LSCI, speckle patterns of consequent frames are correlated in time in the presence of non-moving scatterers, e.g. contains "frozen" speckle patterns. This nonzero speckle pattern correlation is the consequence of the temporal stability of the mutual location of static scatterers, such as skin or tissue, and the laser source. In practice, the result is granular image and false flow rate estimate on static and static scatterer-rich areas. The reason behind it is that static intensity patterns have no proper statistical content to calculate ensemble average. On the other hand, interframe uncorrelated pattern illumination reduces the difference between dynamic and static content, ideally improving the estimation $\propto 1/\sqrt {n}$ regardless of the content [11][24] . In the following, we investigate the use cases of ensemble averaging and how it changes contrast map quality as a function of static and dynamic scatterers ratio.

We are using two premises: the frame by frame illumination change should not diminish the speckle contrast of the individual frames and at the same time, the consecutive patterns should be uncorrelated. Hence, the correlation time $\tau _I$ of the illumination related objective speckle patterns must be significantly larger than the exposure time $T_{exp} \ll \tau _I$, and the decorrelation of the speckle patterns of the consecutive frames must be forced on the whole region of interest.

Different technical solutions could satisfy these conditions: slowly rotating ground glass diffuser [11][25], digital micromirror arrays [26][24], space light phase modulators, tunable laser diodes [27], and other similar techniques that are commonly used in switchable laser speckle reduction solutions [28,29] may be considered for the purpose.

3. Effect of ensemble averaging onstatic and dynamic scatterers

Next, we formulate the combination of intensity speckle patterns of static and dynamic scatterers to model the effect of sample size (number of frames). We adopted the common definition of the fraction of dynamic optical scatterers over total scatterers denoted by $\rho$, and the non-moving scatterers as $(1-\rho )$. The static component is modeled as a single fully developed speckle pattern, the dynamic contribution is described by temporarily integrated distribution [30], and an independent noise source is modeled by Gaussian noise as:

$$I = (1-\rho) S(\lambda_S) + \rho D(k_D,\lambda_D) + W,$$
where $S$ is random variables of exponential distribution with parameter $\lambda _S > 0$ representing static scatterers. Dynamic scatterers are described with the random variables $D$ of a gamma distribution defined by shape-rate parametrization: $k_D,\lambda _D > 0$, where $\lambda _D$ is rate and the $k_D$ is the shape parameter. Camera noise is taken into account as random variable of normal distribution $W \sim \mathcal {N}(\mu _W,\,\sigma _W^{2})$. Any non-zero mean error source such as dark current or readout noise can be modeled by non-zero $\mu _W$ bias and other zero mean sources, like photon or electronic noises, can be described by $\sigma _W^{2}$.

The quality factor is calculated in an alternative way – based on variance and mean. As for the standard deviation of the sample standard deviation no closed-form expression exists. The expression can be approximated by using Taylor series expansion. However, due to the highly nonlinear nature of the problem even the second order approximation introduces significant inaccuracies, the variance based variables are defined by:

$$ K_s^{v} = \frac{\sigma_s'^{2}}{\mu_s'}, \quad \epsilon_K^{v} = \frac{\textrm{Var}{\bigg(}K_s^{v}{\bigg)}}{\textrm{E}{\bigg(}K_s^{v}{\bigg)}}, $$
where $\sigma _s^{2}$ and $\mu _s$ are the sample variance and sample mean, respectively. The $K_s^{v}$ is evaluated similarly to the spatial speckle contrast on a sliding spatial neighborhood of $n_s \times n_s$ pixels than $\epsilon _K^{v}$ is calculated through $n$ frames.

The effect of the number of samples $n$ is used for constructing the variables of the distributions for the two illumination cases. Since variance of sample variance is calculated the dependent variables are extended with $n^{1/4}$. The illumination has no effect on the expected value, thus separate variables are used for the calculation of the variance and the expected value. Using the scaling property of the exponential, the gamma distribution, and the variance, the following change of variables can be made for the different illumination cases.

Stable illumination:

$$\textrm{Var}: \quad \lambda_S' = \frac{\lambda_S}{1-\rho}, \quad \lambda_D' = \frac{\lambda_D}{\rho} n^{1/4}, \quad k_D' = k_D, \mu_W'=\mu_W, \quad \sigma_W' = \frac{\sigma_W}{n^{1/4}},$$
Decorrelated illumination:
$$\textrm{Var}: \quad \lambda_S' = \frac{\lambda_S}{1-\rho}n^{1/4}, \quad \lambda_D' = \frac{\lambda_D}{\rho} n^{1/4},\quad k_D' = k_D \mu_W'=\frac{\mu_W}{n^{1/4}}, \quad \sigma_W' = \frac{\sigma_W}{n^{1/4}}$$
Illumination independent:
$$\textrm{E}: \quad \lambda_S' = \frac{\lambda_S}{1-\rho}, \quad \lambda_S' = \frac{\lambda_D}{\rho},\quad k_D' = k_D \mu_W'=\mu_W, \quad \sigma_W' = \sigma_W,$$
where the $\lambda _S'$, $\lambda _D'$, $k_D'$, $\mu _W'$ and $\sigma _W'$ are the transformed variables. In order to express the quality as a function of sample size, we use moment generating functions [31]:
$$ M_{S}(x) = \frac{\lambda_S'}{\lambda_S' - x}, \quad M_{D}(x) = \Big(\frac{\lambda_D'}{\lambda_D' - x}\Big)^{k_D}, \quad M_{W}(x) = e^{\mu_W x + 0.5(\sigma_W'x)^{2}}, $$
where $M$ is the moment generating function. The k-th moment and the sum of the random variables can be obtained in the following way:
$$ M_{I} = M_{S}(x) M_{D}(x) M_{W}(x) , s$$
$$ E[I^{k}] = \frac{d^{k}}{d x^{k}}M_{I}(x)\Big\rvert_{x=0}. $$
The variance and the mean of the ratio distribution can be approximated using Taylor series expansion [32]. Furthermore, utilizing that the covariance of the sample variance and sample mean is:
$$\textrm{E} \Big(\frac{\sigma_s^{'2}}{\mu_s'}\Big) \approx \frac{\textrm{E}(\sigma_s^{'2})}{\textrm{E}(\mu_s')} - \frac{\textrm{Cov}(\sigma_s^{'2},\mu_s')}{\textrm{E}(\mu_s')^{2}} + \frac{\textrm{Var}(\mu_s')\textrm{E}(\sigma_s^{'2})}{\textrm{E}(\mu_s')^{3}}, $$
$$\textrm{Var}\Big(\frac{\sigma_s^{'2}}{\mu_s'}\Big) \approx \frac{\textrm{Var}(\sigma_s^{'2})}{\textrm{E}(\mu_s')^{2}}+\frac{\textrm{E}(\sigma_s^{'2})^{2}\textrm{Var}(\mu_s')}{\textrm{E}(\mu_s')^{4}}-\frac{2\textrm{Cov}(\sigma_s^{'2},\mu_s')\textrm{E}(\sigma_s^{'2})}{\textrm{E}(\mu_s')^{3}}, $$
$$\textrm{Cov}{\bigg(}\sigma_s^{'2},\mu_s'{\bigg)} = \frac{\mu_3'}{n_s^{2}}, $$
$$\textrm{Var}(\sigma_s^{'2}) = \frac{\mu_4'}{n_s^{2}}-\frac{n_s^{2}-3}{n_s^{2}-1}\frac{\sigma_s^{'4}}{n_s^{2}}, $$
$$\textrm{Var}(\mu_s') = \frac{\sigma_s^{'2}}{n_s^{2}}, $$
$$K_s = \textrm{E}\Big( \frac{\sigma_s}{\mu_s} \Big) \approx \frac{\textrm{E}(\sigma_s')}{\textrm{E}(\mu_s')} - \frac{\textrm{Cov}(\sigma_s',\mu_s')}{\textrm{E}(\mu_s')^{2}} + \frac{\textrm{Var}(\mu_s')\textrm{E}(\sigma_s')}{\textrm{E}(\mu_s')^{3}}, $$
$$\textrm{Cov}{\bigg(}\sigma_s',\mu_s'{\bigg)} = \frac{\mu_3'}{2 \sigma_s' n_s^{2}}, $$
where $\mu _3'$, $\mu _4'$ is the 3rd and 4th central moment and $n_s^{2}$ is the area size. The two most specific limits of the model are given of static and dynamic content only ($\rho =0,1$) for stable illumination, assuming that $\lambda _S$=$\lambda _D$=$\lambda$, the complex relation simplifies to the following:
$$\lim_{\substack{k_D \rightarrow 1 \\ \sigma_W \rightarrow 0 \\ \rho \rightarrow 1 \\ \mu_W \rightarrow 0}} \epsilon_K^{v} = \frac{n_s^{2} -4 n^{1/4}(n_s^{2}-1)+\sqrt{n}(8n_s^{2}-6)-1}{\lambda n^{3/2}(n_s^{2}-1)^{2}},$$
$$\lim_{\substack{k_D \rightarrow 1 \\ \sigma_W \rightarrow 0 \\ \rho \rightarrow 0 \\ \mu_W \rightarrow 0}} \epsilon_K^{v} = \frac{5 n_s^{2} -3}{\lambda (n_s^{2} - 1)^{2}},$$
Equation (18) confirms the expected behavior, that the variance-based quality factor improves with increasing sample number $n$ for a given $n_s^{2}$ area, and shows $\epsilon _K^{v} \propto 1/n$ dependency. Meanwhile, for static content Eq. (19) the quality is sample number $n$ independent at constant illumination, meaning that the averaging does not improve the flatness of the contrast map. On the other hand, the interframe decorrelated illumination makes the samples independent and consequently reduces the variability of the static scatterers as well. In order to validate the model and observe the mixed static, dynamic behavior, a numerical simulation was implemented using the methods and models presented in [33] and [34]. Figure 1 shows a comparison of $\epsilon _K^{v}$ under stable and changing illumination at different $\rho$ values and Fig. 2 demonstrates the effect of the variations of multiple distribution parameters.

 figure: Fig. 1.

Fig. 1. Comparison of the numerical simulations and the model. Quality factor $\epsilon _K^{v}$ is plotted on the left axis, while the spatial contrast $K_s$ with respect to factor $\rho$ on the right axis. The value of parameter $\rho$ defines the fraction of dynamic scatterers and total scatterers ($\rho =1$ only dynamic, $\rho =0$ only static scatterers), $n$ denotes the number of frames acquired, $n_s^{2}$ area was set to constant 400. The numerical simulation, similarly to [33] and [34], was carried out by generating independent random samples for $D$ and $W$ variables, while keeping $S$ variable constant. A sample defines an intensity map as given in Eq. (3). The parameters of the random variables were chosen so that $D$ and $S$ variables had identical mean. The shape parameter was $k_D = 2$, which resulted in a $\sqrt {2}$ spatial contrast in the fully dynamic ($\rho =1$) scenario, the $\mu _W$ and $\sigma _W$ were set to 0.

Download Full Size | PDF

 figure: Fig. 2.

Fig. 2. The model parameters have been varied around a nominal value in order to showcase the effect of the different parameter tendencies on the quality factor and the spatial contrast. The nominal values of the parameters have been based on the validation process described in Sec. 4. The arrows indicate the direction of increase in the parameter values. The camera noise related parameters, namely $\sigma _W$ and $\mu _W$ have opposite effect. Larger mean noise level causes larger spatial contrast degradation, this effect can be attenuated by acquiring several images as the trajectories of the quality factors indicate. Larger $\sigma _W$ values superpose an additional pattern on the image, and increases the spatial contrast. The shape parameter $k_D$ related to the ratio of camera exposure time and decorrelation time, larger values indicate higher ratios, and from a practical point of view a blurred image. Furthermore, the shape parameter doesn’t have an effect on the fully static scenario ($\rho = 0$), as it affects only the distribution of the dynamic variable.

Download Full Size | PDF

The case of temporal contrast calculation must be addressed as well. The temporal contrast calculation based flowmetry is demonstrated to be less sensitive to static scatterers in volume scatterers and flows [17]. Using interframe decorrelated illumination the patterns become statistically independent frame by frame, and as a consequence, the flow rate cannot be estimated with the calculated temporal contrasts. In order to retain the proper temporal statistics and still evaluate the effect of ensemble averaging, we propose a combined experimental method. The temporal contrast is calculated according to Eq. (2) for a short sequence [16] (e.g. 15-20 frames). Then, after providing an independent wavefront phase distribution, another contrast map is gathered. The contrast maps are accumulated in this manner and finally averaged.

4. Demonstrations

Three versatile experiments were conducted to compare the predicted difference between the stable and the ensemble illumination schemes:

  • • Intralipid (Fresenius Kabi AB, Uppsala, Sweden) flow was observed while pumped through a clear channel slide phantom covered with a series of different thickness diffusers.
  • • Blood flow in vessels was visualized in vivo on a surgically prepared mouse brain surface.
  • • The underside of a fresh Ficus benjamina leaf was examined.
The general setup is detailed in Fig. 3. We changed the phase distribution of the scattered field of the laser with a pair of diffusers in which one diffuser was rotated. This arrangement supersede the presented single diffuser arrangements [11][24], in several aspects: i) a doubly scattered speckle system (mm range diameter aperture) approximates better Gaussian statistics, ii) the larger surface structures of one diffuser are suppressed effectively on the illuminated target, iii) the overall speckle pattern count using rotating diffuser, the maximum number of uncorrelated speckle patterns is limited by a single diffuser [24], while the double configuration is still limited, but provides significantly larger independent pattern count [35].

 figure: Fig. 3.

Fig. 3. The optical arrangements of the illumination and imaging system used in the experiments. The examined targets are a channel phantom, a fresh leaf of ficus benjamina, and the surgically prepared mouse brain surface.

Download Full Size | PDF

A laser spot was generated by a collimated and beam expanded single-mode laser diode driven by a constant current source (RLD82PZJ2, 820 nm central wavelength, 220 mW, ROHM Co., Ltd., Japan). The laser diode was placed in a mount with a thermoelectric cooling stage (LDM9T, Thorlabs, Newton, NJ, USA). The images were acquired with a 2x infinity-corrected objective, a 15 cm focal length tube lens, a linear polarizer (Thorlabs, Newton, NJ, USA), and a monochrome camera of 1536x2048 pixels resolution and $3.45x3.45$ µm$^{2}$ pixel size (Basler ACA2040-55um, Basler Vision Technologies, Germany). One fused silica 1500 grit diffuser of the dual pair (DGUV10-1500, Thorlabs, Newton, NJ, USA) was placed in a DC servo motorized rotational stage. The visible area was $4x5$ mm$^{2}$. The speckle size was $s\approx 10.4$ µm. The speed of revolution of the diffuser was set to a slow rotation of $0.05$ rad/s to achieve $\tau _I \approx 50$ ms. We chose exposure times in the range of $2-3$ ms and we set the frame rate to 20 FPS in the demonstrations to satisfy the required premises. The number of emerging temporally independent speckle intensity patterns during the used exposure time is $M = T_{exp}/\tau _I \approx 0.06$, meaning negligible reduction in speckle contrast [35], while each frame is provided an independent pattern. To measure the quality and to quantify the non-uniformity of the speckle contrast maps we use a simpler metric based on Eq. (1):

$$\epsilon_K = \frac{\sigma_K}{\mu_K},$$
where $\epsilon _K$, $\sigma _K$ and $\mu _K$ are the standard deviation and mean values of a region of interest of the calculated contrast [15].

First, we conducted the channel slide experiment to evaluate the effect of ensemble averaging in a controlled environment with known flow speed and changing $\rho$. A syringe pump (SN-50F6, Sino Medical-Device Technology Co., Ltd., China) drove a 10% emulsion of phospholipid stabilized soybean oil (Intralipid) through a µ-Slide I Luer channel slide of 200 µm channel height (ibidi GmbH, Germany) at 2.0 mm/s speed. The slide was placed on a black beam blocking background cavity in order to avoid multiple reflections and provide optical isolation. The entire channel has been covered uniformly with an $\approx 54 \mu m$ thick weak diffuser (Scotch Magic Tape, 3M, USA) with increasing number of layers. During the first measurement the slide was not covered, and one additional layer was applied to the channel slide up step by step after each measurement step, up to 9 layers. In each step, the speckle pattern was recorded for 100-100 frames with an exposure time of $3$ ms with stable and decorrelated illumination as well. The spatial speckle contrast (Eq. (1)) was calculated in a $n_s \times n_s = 7 \times 7$ spatial neighborhood sliding window and then averaged through the frames for $1,2,\ldots ,N$ times. The non-uniformity of the contrast maps is calculated by Eq. (20) using a $n_s \times n_s = 11 \times 11$ spatial window and averaged for each $n$ value. The spatial windows sizes were selected to meet the Nyquist criterion with respect to the speckle size and the spatial contrast pattern. The comparative results are presented in Fig. 4.

 figure: Fig. 4.

Fig. 4. The mean quality factor of the channel slide area is calculated with increasing diffuser coverage (decreasing $\rho$) at constant flow rate 2 mm/s of 10% Intralipid emulsion in a channel slide. Figure (a) shows the quality factor change as a function of the number of averaged contrast maps and $\rho$ with stable illumination. Figure (b) corresponds to the quality factor using decorrelated illumination. Figure (c) maps the mean spatial contrast values of the two illumination cases. In figures (a) and (b) the $\rho = 0$ curves were captured on an ground glass diffuser for comparison.

Download Full Size | PDF

It can be seen in the figure that the increasing static scatterer ratio (decreasing $\rho$) contrast map quality does not improve as the number of frames increases using stable illumination (Fig. 4(a)). On the other hand, as expected, the decorrelated illumination effectively improves the contrast map (Fig. 4(b)). Note, that the data series recorded on static ground glass diffuser of sub speckle sized surface structures and $\rho =0$ follows the $\propto 1/\sqrt {n}$ improvement rate. This is not the case in a structured scatterer: the decreasing slope of improvement is due to the heterogeneous structure of the used diffuser under the quality metric. Comparing the mean spatial contrasts for the two illumination cases Fig. 4(c) suggests that the contrast remains independent of illumination changes. As a conclusion, the expected value of the contrast maps does not change, while it becomes less granulated, limited by the scatterer structure.

The model given in Eqs. (4)–(15) has also been validated on this measurement. The validation is formulated as an identification process with constraints on the $\rho$ factor. Each additional layer of the diffuser is assumed to reduce the $\rho$ factor by a constant $\Delta \rho$ value – as it is illustrated in Fig. 5 –, a fully static $(\rho =0)$ and a fully dynamic ($\rho =1$) sequence is recorded, as well. During the identification beside the model parameters, the $\Delta \rho$ parameter is optimized as follows:

$$\begin{array}{cc} {\begin{array}{c} {\min }\\ {{\Delta }\rho ,{\lambda _D},{\lambda _S},}\\ {{k_D},{\mu _W},{\sigma _W}}\\ {s.t.} \end{array}}&{\frac{{\sum\limits_{n \in \{ 1,100\} } {\sqrt {\frac{1}{{{n_\rho }}}\sum\limits_{k = 1}^{{n_\rho }} ( {\Delta }\mathrm{\epsilon }_K^v(n){)^2}} } }}{{2\mathrm{\epsilon }_{K,max}^v}} + \frac{{\sqrt {\frac{1}{{{n_\rho }}}\sum\limits_{k = 1}^{{n_\rho }} {\Delta } K_s^2} }}{{{K_{s,max}}}}}\\ {}&{\begin{array}{l} {\rho = [0,1 - k {\Delta }\rho ]}\\ {k \in (0,9) \subset {\mathbb Z}}\\ {{k_D} \ge \; 1,} \end{array}} \end{array}$$
where $n_{\rho }=11$ is the number of calculated quality factors and spatial contrasts based on the measurements, $\Delta \epsilon _K^{v}$ denote the difference between the measured and simulated quality factor, while $\Delta K_s$ denote the difference between the spatial contrasts. Furthermore, the quality factors are evaluated at the two extremes ($n=1$ or $n=100$), the contributions from the quality factor and spatial contrast are normalized by their maximum values $\epsilon _{K,max}^{v}$ and $K_{s,max}$, respectively. The results of the identification are summarized in Fig. 5.

 figure: Fig. 5.

Fig. 5. The result of the identification is in good agreement with the measured spatial contrast and quality factor. The identified $\Delta \rho = 0.08$, $\lambda _S = 0.0042$, $\lambda _D = 0.01$, $k_D = 17$, $\mu _W = 0.0015$, $\sigma _W = 0.0016$.

Download Full Size | PDF

Next, we conducted an in vivo experiment to investigate the effect of illumination change. First, we deeply anesthetized a mouse (strain: C57BL6/J) with 25 mg/kg xylazine and 125 mg/kg ketamine in 0.9% NaCl mixture, and then we fixed the head of the animal in a stereotaxic frame (Stoelting Inc, USA). Next we made a 1 cm long sagittal incision on the skull, and cleaned the surface with 3% $H_2O_2$. With a high speed drill we opened the skull and made a 5x5 mm cranial window on the parietal part. The exposure time was $3$ ms and 100-100 frames were recorded with each of the two illumination types. The two speckle contrast maps were calculated using $n_s\times n_s= 7\times 7$ pixels neighborhood and then averaged (Fig. 6).

 figure: Fig. 6.

Fig. 6. Results of the in vivo experiment of the parietal opening of a mouse. Figure (a) brightfield image, (b) averaged spatial speckle contrast map of 100 frames with stable laser illumination, (c) averaged spatial speckle contrast map of the ensemble averaging, (d) cross-sections of (b) and (c) maps along the denoted line. The scale bar is $3\:mm$.

Download Full Size | PDF

After segmenting the brightfield microscopy view to categories of bone/tissue covered parenchyma, uncovered parenchyma, and vessels, we observed significant qualitative improvement ($\sigma _K / \mu _K$ lowered by $45-50\%$) over the cross section at the partially tissue-covered region when decorrelated illumination was used. On the open surface of parenchyma and arteries/large vessels, the decorrelated illumination showed less improvement: $\sigma _K / \mu _K$ is lowered by $15-20\%$ and $10-15\%$, respectively. The later difference can be originated in more heterogeneous and voluminous organization of the parenchyma and deeper light penetration depth compared to the simpler surface vessel structure.

As a representative example for temporal contrast calculation, we examined the midrib and veins of the underside of a fresh Ficus benjamina leaf (Fig. 7(a)). The flow rate varies in the range of a few m/h (<1 mm/s) in the leaves. This slow flow rate does not enable the integral type spatial contrast calculation, thus a moderate sampling rate, short exposure time, and temporal contrast calculation is the appropriate technique to estimate the flow rate [36]. The exposure time was $10$ ms, 100 frames were captured at 20 frames per second and used for temporal contrast calculation. In the reference measurement no diffuser movement was applied as shown in Fig. 7(b). During the modified experiment, the diffuser has been moved after every 20 frame and the contrast calculation was restarted. The result was calculated as the average of the contrasts of the whole sequence (Fig. 7(c)). The quality factor $\sigma _K / \mu _K$ was reduced significantly (Fig. 7(f)), resulting in much less variation. Figure 7(d) cross-section compares a line section of the contrasts near a large and narrow vein. We think that the improvement of this technique in imaging setups, with the limit of ergodicity near static scatterers [37], by partial ensemble averaging is an important observation and worth further investigation. The amount of improvement is similarly to the former experiments, is limited by the volumetric structure of the scatterers.

 figure: Fig. 7.

Fig. 7. Temporal speckle contrast measurement of the underside of a fresh Ficus benjamina leaf. (a) brightfield image, (b) temporal contrast of 100 consecutive frames, (c) average of 5 temporal contrast maps of sequences of 20 frames including diffuser relocation after each sequence. (d) cross-sections of the two (b) and (c) contrast maps. Figures (e-f) show the quality factor of three ROIs for the two cases. The scale bar is $1.5\:mm$.

Download Full Size | PDF

5. Summary

Our conclusions are the followings. The usage of ensemble averaging using decorrelated illumination helps to calculate speckle contrast maps of which local quality becomes independent of the static scatterer content. The technique does not change the relation between broadly revealed speckle contrast and flow speed relations, but significantly reduces the gap between the statistically derived theories and the experienced artifact-prone measurements. We demonstrated that several use cases benefit from this technique, but not all: in the presence of minor static scatterers ($\rho \gtrapprox 0.6$), the ensemble averaging provides diminishing improvements: in biological experiments with direct optical window to the tissues and circulatory system, the reported $\rho$ values varies from $0.75-0.95$ (e.g. [4,22,38]) resulting in an expected $\sigma _K / \mu _K$ contrast map quality factor improvement of $10-15\%$. On the other hand, with tissue covered or embedded veins, arterioles and capillaries, such as skin fold models [3], tooth [17], thin skull covered brain samples [38], and in several agricultural applications the $\rho$ values lie in lower and broader range ($0.15-0.8$) making the ensemble averaging useful for quality enhancement. It provides smoother speckle contrast map and hence better flow rate estimation by $\sigma _K / \mu _K$ improvement of $>20\%$. We found that in structured samples, the quality factor improvement is limited by the static scatterer structure. Notable, that when the deeper tissues are structured heterogeneously, the ensemble averaging also smooths their effect similarly to surface scatterers, which may be advance in blood perfusion imaging. The proposed statistical model is able to substitute numerical simulations, and has been validated with channel slide experiments where it showed good agreement with the measurements.

Funding

Eötvös Loránd Research Network (ELKH KO-40/2020); Nemzeti Kutatási Fejlesztési és Innovációs Hivatal (KDP-1019658).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

References

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

2. P. Li, S. Ni, L. Zhang, S. Zeng, and Q. Luo, “Imaging cerebral blood flow through the intact rat skull with temporal laser speckle imaging,” Opt. Lett. 31(12), 1824–1826 (2006). [CrossRef]  

3. B. Choi, N. M. Kang, and J. Nelson, “Laser speckle imaging for monitoring blood flow dynamics in the in vivo rodent dorsal skin fold model,” Microvasc. Res. 68(2), 143–146 (2004). [CrossRef]  

4. H. Cheng and T. Q. Duong, “Simplified laser-speckle-imaging analysis method and its application to retinal blood flow imaging,” Opt. Lett. 32(15), 2188–2190 (2007). [CrossRef]  

5. P. H. Carvalho, J. B. Barreto, R. A. Braga, and G. F. Rabelo, “Motility parameters assessment of bovine frozen semen by biospeckle laser (bsl) system,” Biosyst. Eng. 102(1), 31–35 (2009). [CrossRef]  

6. R. A. Braga Jr, G. F. Rabelo, L. R. Granato, E. F. Santos, J. C. Machado, R. Arizaga, H. J. Rabal, and M. Trivi, “Detection of fungi in beans by the laser biospeckle technique,” Biosyst. Eng. 91(4), 465–469 (2005). [CrossRef]  

7. G. F. Rabelo, A. M. Enes, R. A. Braga Junior, and I. M. Dal Fabbro, “Frequency response of biospeckle laser images of bean seeds contaminated by fungi,” Biosyst. Eng. 110(3), 297–301 (2011). [CrossRef]  

8. P. Pieczywek, M. Nowacka, M. Dadan, A. Wiktor, K. Rybak, D. Witrowa-Rajchert, and A. Zdunek, “Postharvest monitoring of tomato ripening using the dynamic laser speckle,” Sensors 18(4), 1093 (2018). [CrossRef]  

9. A. Rahmanian, S. A. Mireei, S. Sadri, M. Gholami, and M. Nazeri, “Application of biospeckle laser imaging for early detection of chilling and freezing disorders in orange,” Postharvest Biol. Technol. 162, 111118 (2020). [CrossRef]  

10. M. F. D’Jonsiles, G. E. Galizzi, A. E. Dolinko, M. V. Novas, E. Ceriani Nakamurakare, and C. C. Carmarán, “Optical study of laser biospeckle activity in leaves of jatropha curcas l.: a non-invasive and indirect assessment of foliar endophyte colonization,” Mycol. Prog. 19(4), 339–349 (2020). [CrossRef]  

11. A. C. Völker, P. Zakharov, B. Weber, F. Buck, and F. Scheffold, “Laser speckle imaging with an active noise reduction scheme,” Opt. Express 13(24), 9782–9787 (2005). [CrossRef]  

12. H. Fujii, K. Nohira, Y. Yamamoto, H. Ikawa, and T. Ohura, “Evaluation of blood flow by laser speckle image sensing part 1,” Appl. Opt. 26(24), 5321 (1987). [CrossRef]  

13. R. Arizaga, “Display of local activity using dynamical speckle patterns,” Opt. Eng. 41(2), 287 (2002). [CrossRef]  

14. R. Braga, C. Nobre, A. Costa, T. Sáfadi, and F. da Costa, “Evaluation of activity through dynamic laser speckle using the absolute value of the differences,” Opt. Commun. 284(2), 646–650 (2011). [CrossRef]  

15. J. Qiu, P. Li, W. Luo, J. Wang, H. Zhang, and Q. Luo, “Spatiotemporal laser speckle contrast analysis for blood flow imaging with maximized speckle contrast,” J. Biomed. Opt. 15(1), 016003 (2010). [CrossRef]  

16. 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–565 (2003). [CrossRef]  

17. J. C. Ramirez-San-Juan, C. Regan, B. Coyotl-Ocelotl, and B. Choi, “Spatial versus temporal laser speckle contrast analyses in the presence of static optical scatterers,” J. Biomed. Opt. 19(10), 106009 (2014). [CrossRef]  

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

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

20. T. Dragojević, D. Bronzi, H. M. Varma, C. P. Valdes, C. Castellvi, F. Villa, A. Tosi, C. Justicia, F. Zappa, and T. Durduran, “High-speed multi-exposure laser speckle contrast imaging with a single-photon counting camera,” Biomed. Opt. Express 6(8), 2865–2876 (2015). [CrossRef]  

21. C. Wang, Z. Cao, X. Jin, W. Lin, Y. Zheng, B. Zeng, and M. Xu, “Robust quantitative single-exposure laser speckle imaging with true flow speckle contrast in the temporal and spatial domains,” Biomed. Opt. Express 10(8), 4097–4114 (2019). [CrossRef]  

22. D. D. Postnov, J. Tang, S. E. Erdener, K. Kiliç, and D. A. Boas, “Dynamic light scattering imaging,” Sci. Adv. 6(45), eabc4628 (2020). [CrossRef]  

23. M. Siket, I. Jánoki, K. Demeter, M. Szabó, and P. Földesy, “Time varied illumination laser speckle contrast imaging,” Opt. Lett. 46(4), 713–716 (2021). [CrossRef]  

24. L. Schweickhardt, A. Tausendfreund, D. Stöbener, and A. Fischer, “Noise reduction in high-resolution speckle displacement measurements through ensemble averaging,” Appl. Opt. 60(7), 1871–1880 (2021). [CrossRef]  

25. T. Stangner, H. Zhang, T. Dahlberg, K. Wiklund, and M. Andersson, “Step-by-step guide to reduce spatial coherence of laser light using a rotating ground glass diffuser,” Appl. Opt. 56(19), 5427–5435 (2017). [CrossRef]  

26. M. N. Akram, Z. Tong, G. Ouyang, X. Chen, and V. Kartashov, “Laser speckle reduction due to spatial and angular diversity introduced by fast scanning micromirror,” Appl. Opt. 49(17), 3297–3304 (2010). [CrossRef]  

27. H. Nasim and Y. Jamil, “Recent advancements in spectroscopy using tunable diode lasers,” Laser Phys. Lett. 10(4), 043001 (2013). [CrossRef]  

28. X. Chen, Ø. Svensen, and M. N. Akram, “Speckle reduction in laser projection using a dynamic deformable mirror,” Opt. Express 22(9), 11152–11166 (2014). [CrossRef]  

29. M. Elbaum, M. Greenebaum, and M. King, “A wavelength diversity technique for reduction of speckle size,” Opt. Commun. 5(3), 171–174 (1972). [CrossRef]  

30. J. W. Goodman, “Statistical properties of laser speckle patterns,” in Laser speckle and related phenomena, (Springer, 1975), pp. 9–75.

31. K. Ramachandran and C. Tsokos, Mathematical Statistics with Applications (Elsevier Science, 2009).

32. J. Rice, Mathematical Statistics and Data Analysis, Advanced series (Cengage Learning, 2006).

33. L. Song, Z. Zhou, X. Wang, X. Zhao, and D. S. Elson, “Simulation of speckle patterns with pre-defined correlation distributions,” Biomed. Opt. Express 7(3), 798–809 (2016). [CrossRef]  

34. H. Spahr, C. Pfäffle, S. Burhan, L. Kutzner, F. Hilge, G. Hüttmann, and D. Hillmann, “Phase-sensitive interferometry of decorrelated speckle patterns,” Sci. Rep. 9(1), 11748 (2019). [CrossRef]  

35. D. Li, D. P. Kelly, and J. T. Sheridan, “Speckle suppression by doubly scattering systems,” Appl. Opt. 52(35), 8617–8626 (2013). [CrossRef]  

36. P. M. Pieczywek, J. Cybulska, A. Zdunek, and A. Kurenda, “Exponentially smoothed fujii index for online imaging of biospeckle spatial activity,” Comput. Electron. Agric. 142, 70–78 (2017). [CrossRef]  

37. P. Zakharov, “Ergodic and non-ergodic regimes in temporal laser speckle imaging,” Opt. Lett. 42(12), 2299–2301 (2017). [CrossRef]  

38. A. B. Parthasarathy, S. S. Kazmi, and A. K. Dunn, “Quantitative imaging of ischemic stroke through thinned skull in mice with multi exposure speckle imaging,” Biomed. Opt. Express 1(1), 246–259 (2010). [CrossRef]  

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (7)

Fig. 1.
Fig. 1. Comparison of the numerical simulations and the model. Quality factor $\epsilon _K^{v}$ is plotted on the left axis, while the spatial contrast $K_s$ with respect to factor $\rho$ on the right axis. The value of parameter $\rho$ defines the fraction of dynamic scatterers and total scatterers ($\rho =1$ only dynamic, $\rho =0$ only static scatterers), $n$ denotes the number of frames acquired, $n_s^{2}$ area was set to constant 400. The numerical simulation, similarly to [33] and [34], was carried out by generating independent random samples for $D$ and $W$ variables, while keeping $S$ variable constant. A sample defines an intensity map as given in Eq. (3). The parameters of the random variables were chosen so that $D$ and $S$ variables had identical mean. The shape parameter was $k_D = 2$, which resulted in a $\sqrt {2}$ spatial contrast in the fully dynamic ($\rho =1$) scenario, the $\mu _W$ and $\sigma _W$ were set to 0.
Fig. 2.
Fig. 2. The model parameters have been varied around a nominal value in order to showcase the effect of the different parameter tendencies on the quality factor and the spatial contrast. The nominal values of the parameters have been based on the validation process described in Sec. 4. The arrows indicate the direction of increase in the parameter values. The camera noise related parameters, namely $\sigma _W$ and $\mu _W$ have opposite effect. Larger mean noise level causes larger spatial contrast degradation, this effect can be attenuated by acquiring several images as the trajectories of the quality factors indicate. Larger $\sigma _W$ values superpose an additional pattern on the image, and increases the spatial contrast. The shape parameter $k_D$ related to the ratio of camera exposure time and decorrelation time, larger values indicate higher ratios, and from a practical point of view a blurred image. Furthermore, the shape parameter doesn’t have an effect on the fully static scenario ($\rho = 0$), as it affects only the distribution of the dynamic variable.
Fig. 3.
Fig. 3. The optical arrangements of the illumination and imaging system used in the experiments. The examined targets are a channel phantom, a fresh leaf of ficus benjamina, and the surgically prepared mouse brain surface.
Fig. 4.
Fig. 4. The mean quality factor of the channel slide area is calculated with increasing diffuser coverage (decreasing $\rho$) at constant flow rate 2 mm/s of 10% Intralipid emulsion in a channel slide. Figure (a) shows the quality factor change as a function of the number of averaged contrast maps and $\rho$ with stable illumination. Figure (b) corresponds to the quality factor using decorrelated illumination. Figure (c) maps the mean spatial contrast values of the two illumination cases. In figures (a) and (b) the $\rho = 0$ curves were captured on an ground glass diffuser for comparison.
Fig. 5.
Fig. 5. The result of the identification is in good agreement with the measured spatial contrast and quality factor. The identified $\Delta \rho = 0.08$, $\lambda _S = 0.0042$, $\lambda _D = 0.01$, $k_D = 17$, $\mu _W = 0.0015$, $\sigma _W = 0.0016$.
Fig. 6.
Fig. 6. Results of the in vivo experiment of the parietal opening of a mouse. Figure (a) brightfield image, (b) averaged spatial speckle contrast map of 100 frames with stable laser illumination, (c) averaged spatial speckle contrast map of the ensemble averaging, (d) cross-sections of (b) and (c) maps along the denoted line. The scale bar is $3\:mm$.
Fig. 7.
Fig. 7. Temporal speckle contrast measurement of the underside of a fresh Ficus benjamina leaf. (a) brightfield image, (b) temporal contrast of 100 consecutive frames, (c) average of 5 temporal contrast maps of sequences of 20 frames including diffuser relocation after each sequence. (d) cross-sections of the two (b) and (c) contrast maps. Figures (e-f) show the quality factor of three ROIs for the two cases. The scale bar is $1.5\:mm$.

Equations (21)

Equations on this page are rendered with MathJax. Learn more.

K s = σ ( I s ) I s ,
K t = σ ( I t ) I t = I t 2 I t 2 I t 2 ,
I = ( 1 ρ ) S ( λ S ) + ρ D ( k D , λ D ) + W ,
K s v = σ s 2 μ s , ϵ K v = Var ( K s v ) E ( K s v ) ,
Var : λ S = λ S 1 ρ , λ D = λ D ρ n 1 / 4 , k D = k D , μ W = μ W , σ W = σ W n 1 / 4 ,
Var : λ S = λ S 1 ρ n 1 / 4 , λ D = λ D ρ n 1 / 4 , k D = k D μ W = μ W n 1 / 4 , σ W = σ W n 1 / 4
E : λ S = λ S 1 ρ , λ S = λ D ρ , k D = k D μ W = μ W , σ W = σ W ,
M S ( x ) = λ S λ S x , M D ( x ) = ( λ D λ D x ) k D , M W ( x ) = e μ W x + 0.5 ( σ W x ) 2 ,
M I = M S ( x ) M D ( x ) M W ( x ) , s
E [ I k ] = d k d x k M I ( x ) | x = 0 .
E ( σ s 2 μ s ) E ( σ s 2 ) E ( μ s ) Cov ( σ s 2 , μ s ) E ( μ s ) 2 + Var ( μ s ) E ( σ s 2 ) E ( μ s ) 3 ,
Var ( σ s 2 μ s ) Var ( σ s 2 ) E ( μ s ) 2 + E ( σ s 2 ) 2 Var ( μ s ) E ( μ s ) 4 2 Cov ( σ s 2 , μ s ) E ( σ s 2 ) E ( μ s ) 3 ,
Cov ( σ s 2 , μ s ) = μ 3 n s 2 ,
Var ( σ s 2 ) = μ 4 n s 2 n s 2 3 n s 2 1 σ s 4 n s 2 ,
Var ( μ s ) = σ s 2 n s 2 ,
K s = E ( σ s μ s ) E ( σ s ) E ( μ s ) Cov ( σ s , μ s ) E ( μ s ) 2 + Var ( μ s ) E ( σ s ) E ( μ s ) 3 ,
Cov ( σ s , μ s ) = μ 3 2 σ s n s 2 ,
lim k D 1 σ W 0 ρ 1 μ W 0 ϵ K v = n s 2 4 n 1 / 4 ( n s 2 1 ) + n ( 8 n s 2 6 ) 1 λ n 3 / 2 ( n s 2 1 ) 2 ,
lim k D 1 σ W 0 ρ 0 μ W 0 ϵ K v = 5 n s 2 3 λ ( n s 2 1 ) 2 ,
ϵ K = σ K μ K ,
min Δ ρ , λ D , λ S , k D , μ W , σ W s . t . n { 1 , 100 } 1 n ρ k = 1 n ρ ( Δ ϵ K v ( n ) ) 2 2 ϵ K , m a x v + 1 n ρ k = 1 n ρ Δ K s 2 K s , m a x ρ = [ 0 , 1 k Δ ρ ] k ( 0 , 9 ) Z k D 1 ,
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.