1.

Introduction

A speckle pattern is created when coherent light is reflected from a rough surface or transmitted through a scattering media. This pattern is made up of regions with high and low intensity and depends on the distribution of scatterers. Static scatterers create a static speckle pattern, while dynamic scatterers create a dynamic speckle pattern. Speckle patterns appear in images collected with many detection systems, including radar,^1,2 microscopy,^3–5 astronomy,^6–8 and ultrasound.^9,10 Biological samples produce speckle patterns with different dynamics in different regions, such as those associated with blood vessels. Regions rich in blood vessels, like skin, are of particular interest to monitor skin lesions, such as wounds, ulcers, and tumors. It has also been used to assess the efficacy of various treatments, such as topical creams and laser therapy. Additionally, speckle pattern imaging can be used to study the effects of aging and diseases, such as diabetes, on skin microcirculation. When recorded by a camera with finite exposure time, dynamic speckle patterns are blurred, and this local blurring is related to the movement of scatterers, such as blood flow, and can be measured as speckle contrast ($K$).

Speckle contrast analysis is the basis of laser speckle contrast imaging (LSCI), a low-cost, noninvasive technique used in different areas of medicine, such as dentistry,^11,12 ophthalmology,^13,14 dermatology,^15,16 neurobiology,^17,18 and neuroscience.^19,20 Common uses of LSCI are to measure blood flow^21,22 and perform intraoperative assessment of tissue perfusion.²³

For LSCI, the propagation of coherent light through scattering media produces a speckle pattern captured with a digital camera. If the medium contains moving scatterers, the speckle pattern fluctuates proportionally to the scatterers’ speed. If the exposure time is longer than the correlation time of the backscattered light, the speckle visibility will be reduced. In 1981, Fercher and Briers first measured blood flow speed using LSCI²⁴ and derived the following expression for the contrast $K$:

Since this seminal paper, the LSCI technique has evolved and therefore the contrast calculation too. For example, Lemieux and Durian²⁵ added a correction factor $\beta$ (which depends on the ratio of the pixel and speckle grain size) to the Siegert relation that links the intensity dynamics associated with scatterer motion with the electric field dynamics,²⁶ resulting in $g_2(\tau )=1+\beta {|g_1(\tau )|}^2$, where $g_1(\tau )$ is the electric field autocorrelation function and $g_2(\tau )$ is the intensity autocorrelation function. Bandyopadhyay et al.²⁷ corrected an over simplification on the derivation of Eq. (1), applying a double integral of $\tau =t_1-t_2$, obtaining

Bandyopadhyay et al.²⁷ proposed a technique to measure $\beta$ experimentally. Accounting for the finite size of the camera pixels²⁹ led to an analytical expression for $\beta$, assuming Gaussian illumination:

With the derivation of the equations above, the correlation between neighboring pixels is neglected, and hence are valid only for small speckle area compared to the pixel area ($M\ge 1$). When the speckle area is larger than the pixel area ($M<1$), these models do not fit experimental data well. Here we present a mathematical model that considers correlations between neighboring pixels to improve the agreement between theoretical and experimental contrast measurements for any size of speckle and pixel areas $(M>0)$.

2.

Theory

2.1.

Spatial Correlation Model

Several computational algorithms^30,31 exist to calculate the speckle contrast from a single raw speckle image or set of images. Application of spatial statistics^24,32 to a single image is probably the most common approach. This analysis typically implements a sliding square window to calculate the local speckle contrast as the standard deviation of pixel intensities divided by the mean intensity. However, the process of the sliding window is not taken into in the current theoretical model (any of the above mentioned models), which can lead to an inaccurate contrast calculation, and therefore to an inaccurate blood flow calculation. Indeed, it is necessary to consider not only the contribution of the physical size of the camera’s pixel (autocorrelation) but also the effect of the size of the sliding window to improve the accuracy of $\beta$ and therefore of the contrast.

Skipetrov et al.³³ employed a sliding window of $3\times 3$ pixels. In this work, we expand on their analysis to extend it to a sliding window of any size to measure an accurate blood flow speed. We show that the mathematical model considering the statistical analysis of the sliding window agrees well with experimental data.

We assume a speckle pattern $I(\overrightarrow{r})$ recorded by a digital camera, where $\overrightarrow{r}=x\hat{i}+y\widehat{j}$ is the position vector on the CCD image. Here we consider a sliding window of $\sqrt{N}\times \sqrt{N}$ pixels, with $N$ being an odd number. We define the speckle spot size as comprising a subregion of $(2p+1)\times (2p+1)$ pixels in size, where $p\in {\mathbb{Z}}^{+}$ and $2p+1\le \sqrt{N}$. Within this subregion, the pixels are assumed to be spatially correlated. Figure 1(a) shows a schematic representation of the sliding window (pink rectangles) and one possible example of a correlation subregion (blue triangles).

Fig. 1

(a) Schematic representation of the sliding window of ($\sqrt{N}\times \sqrt{N}$) pixels (pink) and the correlation subregion of $(2p+1{)}^2$ (blue). In this example, the sliding window of size $\sqrt{N}=11$, and correlation subregion of size $2p+1=7$ is shown. (b) Schematic representation of the correlation region of $(2p+1{)}^2$ pixels, with a close-up view of the first quadrant, and a chess board showing the movement of a knight. Some pixels are labeled with their coordinates, and they are colored according to their type (see text for details).

Figure 1(b) shows the corresponding coordinates $(\eta ,\xi )$ for each pixel of the sliding window, where $\eta ,\xi \in \{0,\pm 1,\pm 2\dots ,\pm \frac{\sqrt{N}-1}{2}\}$ and a correlation subregion containing different types of correlations: central, lateral, diagonal, and knight. It is important to mention that all pixels within the subregion are correlated and all of them are considered on calculations. However, due to symmetry, many of the correlations are repeated several times. For example, since the lateral pixels $(0,\eta )$, $(\eta ,0)$, $(0,-\eta )$, and $(-\eta ,0)$ are all at the same distance from the central pixel of the sliding window, they are equally correlated with the central pixel, and therefore, this correlation value is repeated four times.

Considering a pixel with area $a^2$, the analytic expression for the speckle contrast calculated over an arbitrary sliding window size of $\sqrt{N}\times \sqrt{N}$ pixels, and hence an arbitrary correlation subregion of $(2p+1)\times (2p+1)$ pixels may be written as (details can be found in Appendix):

Eq. (6)

$$\begin{gathered} K_s^2(N,p)=\frac{1}{{\mu }_{\mathrm{0,0}}}-\frac{1}{N(N-1)}\left(4\sum _{\eta =1}^p\right(\sqrt{N}-\eta )\sqrt{N}\frac{1}{{\mu }_{\eta ,0}}+4\sum _{\eta =1}^p{(\sqrt{N}-\eta )}^2\frac{1}{{\mu }_{\eta ,\eta }} \\ +8\sum _{\xi =1}^{p-1}\sum _{\eta =\xi +1}^p(\sqrt{N}-\eta \left)\right(\sqrt{N}-\xi \left)\frac{1}{{\mu }_{\eta ,\xi }}\right), \end{gathered}$$

where each summation inside the parentheses stands for a diagonal, lateral, or knight-type correlation, respectively. As described in Appendix, Eq. (6) is valid for $p\ge 2$, i.e., a $5\times 5$ subregion. For $p=1$, there is no knight correlation, and thus the double summation does not exist. For $p=0$, there is no lateral, diagonal, or knight-type correlation, and thus the terms in parentheses equal zero. Here we consider only $p\ge 2$. $1/{\mu }_{\eta ,\xi }$ is the correlation factor of the pixel with coordinates $(\eta ,\xi )$ and with the central pixel $(\mathrm{0,0})$ of the sliding window, given by

Eq. (7)

$$\frac{1}{{\mu }_{\eta ,\xi }}=\frac{1}{a^4}{\int }_{a(2\xi -1)/2}^{a(2\xi +1)/2}{\int }_{a(2\eta -1)/2}^{a(2\eta +1)/2}{\int }_{-a/2}^{a/2}{\int }_{-a/2}^{a/2}{g}_1^2(x-x^{\prime },y-y^{\prime })\mathrm{d}x\mathrm{d}y\mathrm{d}{x}^{\prime }\mathrm{d}{y}^{\prime },$$

where $g_1$ is a first-order correlation function. At this point, we have not considered the correlation function $(g_1^2(\overrightarrow{r}-{\overrightarrow{r}}^{\prime }))$; this will be the topic of the next section. The factor $\beta$ in Eq. (4) corresponds to the spatial contrast $K_s^2(N,p)$ (i.e., $\beta =K_s^2(N,p)$) for any window size, while Eq. (5), $K_s^2(N,p)$ takes into account not only the finite dimension of the pixel but its possible spatial correlation.

Table 1 shows expressions for the spatial contrast using Eq. (6) for the first four different dimensions of the sliding window.

Table 1

Expressions for Eq. (6) for sliding window sizes of 1×1 (i.e., autocorrelation), 3×3, 5×5, and 7×7  pixels.

$K_s^2(1^2,0)$	$\frac{1}{{\mu }_{\mathrm{0,0}}}$
$K_{\mathrm{s}}^2(3^2,1)$	$\frac{1}{{\mu }_{\mathrm{0,0}}}+(-\frac{1}{3}\frac{1}{{\mu }_{\mathrm{1,0}}})+(-\frac{2}{9}\frac{1}{{\mu }_{\mathrm{1,1}}})$
$K_{\mathrm{s}}^2(5^2,2)$	$\frac{1}{{\mu }_{\mathrm{0,0}}}+(-\frac{2}{15}\frac{1}{{\mu }_{\mathrm{1,0}}}-\frac{1}{10}\frac{1}{{\mu }_{\mathrm{2,0}}})+(-\frac{8}{75}\frac{1}{{\mu }_{\mathrm{1,1}}}-\frac{3}{50}\frac{1}{{\mu }_{\mathrm{2,2}}})+(-\frac{4}{25}\frac{1}{{\mu }_{\mathrm{2,1}}})$
$K_{\mathrm{s}}^2(7^2,3)$	$\frac{1}{{\mu }_{\mathrm{0,0}}}+(-\frac{1}{14}\frac{1}{{\mu }_{\mathrm{1,0}}}-\frac{5}{84}\frac{1}{{\mu }_{\mathrm{2,0}}}-\frac{1}{21}\frac{1}{{\mu }_{\mathrm{3,0}}})+(-\frac{3}{49}\frac{1}{{\mu }_{\mathrm{1,1}}}-\frac{25}{588}\frac{1}{{\mu }_{\mathrm{2,2}}}-\frac{4}{147}\frac{1}{{\mu }_{\mathrm{3,3}}})+(-\frac{5}{49}\frac{1}{{\mu }_{\mathrm{2,1}}}-\frac{4}{49}\frac{1}{{\mu }_{\mathrm{3,1}}}-\frac{10}{147}\frac{1}{{\mu }_{\mathrm{3,2}}})$

As it seen in Table 1, the contribution of each lateral, diagonal, and knight correlation decreases with increasing subregion size.

2.2.

Gaussian Velocity Distribution Function

For a Gaussian-shaped laser beam, the spatial correlation function also is Gaussian in shape.

2.2.1.

Analytical solutions

The Gaussian correlation function $g_1^2(\overrightarrow{r}-{\overrightarrow{r}}^{\prime })=\exp [-\frac{\pi }{2}\frac{{(x-x^{\prime })}^2+{(y-y^{\prime })}^2}{A_c}]$^33,34 may be written using a change of variables $w=x-x^{\prime }$, $\overline{w}=(x+x^{\prime })/2$, $u=y-y^{\prime }$, and $\overline{u}=(y+y^{\prime })/2$:

Eq. (8)

$$g_1(w,u)=\exp [-\frac{\pi }{2}\frac{(w^2+u^2)}{A_c}],$$

where $A_c$ is the correlation area (i.e., speckle size). For a speckle area with radius $b$:

Eq. (9)

$$g_1^2(w,u)={\left(\exp \right[-\frac{\pi }{2}\frac{(w^2+u^2)}{\pi b^2}\left]\right)}^2=\exp [-(w^2+u^2)/b^2].$$

Substituting Eq. (9) into Eq. (7), we obtain

Eq. (10)

$$\frac{1}{{\mu }_{\eta ,\xi }}=\frac{1}{4{\pi }^2{M}^2}\times (e^{-\pi M(\eta -1{)}^2}-2e^{-\pi M{\eta }^2}+e^{-\pi M(\eta +1{)}^2}+\pi \sqrt{M}((\eta -1)\mathrm{erf}[\sqrt{\pi M}(\eta -1)]-2\eta \mathrm{erf}\left[\sqrt{\pi M}\eta \right]+(\eta +1)\mathrm{erf}[\sqrt{\pi M}(\eta +1)]\left)\right)\times (e^{-\pi M{(\xi -1)}^2}-2e^{-\pi M{\xi }^2}+e^{-\pi M{(\xi +1)}^2}+\pi \sqrt{M}((\xi -1)\mathrm{erf}[\sqrt{\pi M}(\xi -1)]-2\xi \mathrm{erf}\left(\sqrt{\pi M}\xi \right)+(\xi +1)\mathrm{erf}[\sqrt{\pi M}(\xi +1)]\left)\right),$$

where $M=a^2/\pi b^2$ is the ratio of the pixel area and the speckle area. Equation (10) provides an expression for speckle contrast that is similar to that described in section 4.6 of Ref. 35.

Equations (6) and (10) are programmed and available in MATLAB and Mathematica (see Code and Data Availability).

One particular case is when $p=0$ (i.e., autocorrelation), corresponding to a spatial correlation subregion of only a single pixel. For this case, we obtain the well-established analytic expression for spatial speckle contrast:²⁹

Eq. (11)

$$K_s^2(N,p=0)=\frac{1}{{\mu }_{\mathrm{0,0}}}=\frac{1}{{\pi }^2{M}^2}{(-1+e^{-\pi M}+\pi \sqrt{M}\mathrm{erf}[\sqrt{\pi M}\left]\right)}^2.$$

Equation (11) is identical to the correction factor $\beta$ [see Eq. (5)] proposed in Refs. 19 and 25.

Table 2 shows the first $1/{\mu }_{\eta ,\xi }$ factors needed to calculate the speckle contrast for $1\times 1$, $3\times 3$, $5\times 5$, and $7\times 7$ sliding windows.

Table 2

Correlation factors calculated using Eq. (10) needed to calculate the speckle contrast of a 7×7 sliding window.

Central
$\frac{1}{{\mu }_{\mathrm{0,0}}}$	$\frac{1}{M^2{\pi }^2}{(-1+{\mathrm{e}}^{-M\pi }+\sqrt{M}\pi \mathrm{erf}(\sqrt{M}\sqrt{\pi }\left)\right)}^2$
Lateral
$\frac{1}{{\mu }_{\mathrm{1,0}}}$	$$\begin{gathered} \frac{1}{2M^2{\pi }^2}(-1+{\mathrm{e}}^{-M\pi }+\sqrt{M}\pi \mathrm{erf}(\sqrt{M}\sqrt{\pi }\left)\right)(1+{\mathrm{e}}^{-4M\pi }-2{\mathrm{e}}^{-M\pi } \\ -2\sqrt{M}\pi (\mathrm{erf}\left(\sqrt{M}\sqrt{\pi }\right)-\mathrm{erf}\left(2\sqrt{M}\sqrt{\pi }\right)\left)\right) \end{gathered}$$
$\frac{1}{{\mu }_{\mathrm{2,0}}}$	$$\begin{gathered} \frac{1}{2M^2{\pi }^2}(-1+{\mathrm{e}}^{-M\pi }+\sqrt{M}\pi \mathrm{erf}(\sqrt{M}\sqrt{\pi }\left)\right)({\mathrm{e}}^{-9M\pi }-2{\mathrm{e}}^{-4M\pi }+{\mathrm{e}}^{-M\pi } \\ +\sqrt{M}\pi (\mathrm{erf}\left(\sqrt{M}\sqrt{\pi }\right)-4\mathrm{erf}\left(2\sqrt{M}\sqrt{\pi }\right)+3\mathrm{erf}\left(3\sqrt{M}\sqrt{\pi }\right)\left)\right) \end{gathered}$$
$\frac{1}{{\mu }_{\mathrm{3,0}}}$	$$\begin{gathered} \frac{1}{2M^2{\pi }^2}(-1+{\mathrm{e}}^{-M\pi }+\sqrt{M}\pi \mathrm{erf}(\sqrt{M}\sqrt{\pi }\left)\right)({\mathrm{e}}^{-16M\pi }-2{\mathrm{e}}^{-9M\pi }+{\mathrm{e}}^{-4M\pi } \\ +2\sqrt{M}\pi (\mathrm{erf}\left(2\sqrt{M}\sqrt{\pi }\right)-3\mathrm{erf}\left(3\sqrt{M}\sqrt{\pi }\right)+2\mathrm{erf}\left(4\sqrt{M}\sqrt{\pi }\right)\left)\right) \end{gathered}$$
Diagonal
$\frac{1}{{\mu }_{\mathrm{1,1}}}$	$\frac{1}{4M^2{\pi }^2}(1+{\mathrm{e}}^{-4M\pi }-2{\mathrm{e}}^{-M\pi }-2\sqrt{M}\pi {\left(\mathrm{erf}\right(\sqrt{M}\sqrt{\pi })-\mathrm{erf}(2\sqrt{M}\sqrt{\pi }\left)\right))}^2$
$\frac{1}{{\mu }_{\mathrm{2,2}}}$	$\frac{1}{4M^2{\pi }^2}({\mathrm{e}}^{-9M\pi }-2{\mathrm{e}}^{-4M\pi }+{\mathrm{e}}^{-M\pi }+\sqrt{M}\pi {\left(\mathrm{erf}\right(\sqrt{M}\sqrt{\pi })-4\mathrm{erf}(2\sqrt{M}\sqrt{\pi })+3\mathrm{erf}(3\sqrt{M}\sqrt{\pi }\left)\right))}^2$
$\frac{1}{{\mu }_{\mathrm{3,3}}}$	$\frac{1}{4M^2{\pi }^2}({\mathrm{e}}^{-16M\pi }-2{\mathrm{e}}^{-9M\pi }+{\mathrm{e}}^{-4M\pi }+2\sqrt{M}\pi (\mathrm{erf}\left(2\sqrt{M}\sqrt{\pi }\right)-3\mathrm{erf}\left(3\sqrt{M}\sqrt{\pi }\right)+2\mathrm{erf}\left(4\sqrt{M}\sqrt{\pi }\right)){)}^2$
Knight
$\frac{1}{{\mu }_{\mathrm{2,1}}}$	$$\begin{gathered} \frac{1}{4M^2{\pi }^2}(1+{\mathrm{e}}^{-4M\pi }-2{\mathrm{e}}^{-M\pi }-2\sqrt{M}\pi (\mathrm{erf}\left(\sqrt{M}\sqrt{\pi }\right)-\mathrm{erf}\left(2\sqrt{M}\sqrt{\pi }\right)\left)\right)({\mathrm{e}}^{-9M\pi }-2{\mathrm{e}}^{-4M\pi }+{\mathrm{e}}^{-M\pi } \\ +\sqrt{M}\pi (\mathrm{erf}\left(\sqrt{M}\sqrt{\pi }\right)-4\mathrm{erf}\left(2\sqrt{M}\sqrt{\pi }\right)+3\mathrm{erf}\left(3\sqrt{M}\sqrt{\pi }\right)\left)\right) \end{gathered}$$
$\frac{1}{{\mu }_{\mathrm{3,1}}}$	$$\begin{gathered} \frac{1}{4M^2{\pi }^2}(1+{\mathrm{e}}^{-4M\pi }-2{\mathrm{e}}^{-M\pi }-2\sqrt{M}\pi (\mathrm{erf}\left(\sqrt{M}\sqrt{\pi }\right)-\mathrm{erf}\left(2\sqrt{M}\sqrt{\pi }\right)\left)\right)({\mathrm{e}}^{-16M\pi }-2{\mathrm{e}}^{-9M\pi }+{\mathrm{e}}^{-4M\pi } \\ +2\sqrt{M}\pi (\mathrm{erf}\left(2\sqrt{M}\sqrt{\pi }\right)-3\mathrm{erf}\left(3\sqrt{M}\sqrt{\pi }\right)+2\mathrm{erf}\left(4\sqrt{M}\sqrt{\pi }\right)\left)\right) \end{gathered}$$
$\frac{1}{{\mu }_{\mathrm{3,2}}}$	$$\begin{gathered} \frac{1}{4M^2{\pi }^2}({\mathrm{e}}^{-9M\pi }-2{\mathrm{e}}^{-4M\pi }+{\mathrm{e}}^{-M\pi }+\sqrt{M}\pi (\mathrm{erf}\left(\sqrt{M}\sqrt{\pi }\right)-4\mathrm{erf}\left(2\sqrt{M}\sqrt{\pi }\right) \\ +3\mathrm{erf}\left(3\sqrt{M}\sqrt{\pi }\right)\left)\right)({\mathrm{e}}^{-16M\pi }-2{\mathrm{e}}^{-9M\pi }+{\mathrm{e}}^{-4M\pi } \\ +2\sqrt{M}\pi (\mathrm{erf}\left(2\sqrt{M}\sqrt{\pi }\right)-3\mathrm{erf}\left(3\sqrt{M}\sqrt{\pi }\right)+2\mathrm{erf}\left(4\sqrt{M}\sqrt{\pi }\right)\left)\right) \end{gathered}$$

2.2.2.

Numerical solution

We integrated Eq. (7) numerically using Mathematica 11 and compared the resultant solution to the analytical solution in Eq. (10). A good agreement between the analytic and numerical solution is observed (Fig. 2).

Fig. 2

Comparison between the analytical expression (solid lines) using Eq. (10) and a numerical solution (symbols) of Eq. (7) using Mathematica 11.

The parameters $1/{\mu }_{\mathrm{0,0}}$, $1/{\mu }_{\mathrm{1,0}}$, and $1/{\mu }_{\mathrm{1,1}}$ are sufficient to calculate the contrast $K$ for a spatial correlation subregion of $3\times 3\text{pixels}$.³⁶ For a subregion of $5\times 5\text{pixels}$, we also should include the terms $1/{\mu }_{\mathrm{2,0}}$, $1/{\mu }_{\mathrm{2,1}}$, and $1/{\mu }_{\mathrm{2,2}}$.³⁷ Now it is possible to calculate the speckle contrast using Eq. (10) for any arbitrary odd size of sliding window ($\sqrt{N}$) and subregion ($2p+1$).

Figure 3 shows the plots of the analytical (solid lines) and numerical (symbols) contrast for a correlation subregion of $1\times 1$ (i.e., no spatial correlation), $3\times 3$, and $5\times 5$ pixels using Eq. (6). When the speckle size is larger than the pixel size (i.e., $M<1$), each speckle grain spans multiple pixels, resulting in spatial correlation among adjacent pixels within the sliding window. Due to this correlation, the standard deviation $\sigma$ of intensity values within the sliding window decreases for $M<1$, resulting in an associated decrease in $K$. When the speckle size is smaller than the pixel size (i.e., $M>1$), each pixel contains intensity contributions from multiple speckle grains, leading to a lack of spatial correlation among adjacent pixels within the sliding window. In addition, the integration of multiple grains per pixel leads to a reduction in intensity variance in the speckle pattern.

Fig. 3

Solid lines show the contrast calculated using Eq. (6) for a correlation subregion of $1\times 1$ (no spatial correlation), $3\times 3$, and $5\times 5\text{pixels}$, in blue, pink, and black, respectively. The circle, triangle, and square marks show the numerical solution using Mathematica 11 for spatial correlation submatrices of $1\times 1$, $3\times 3$, and $5\times 5\text{pixels}$. Good agreement is observed between the numerical and analytical expressions.

In Fig. 3, it is show that for each size of sliding window, there is a different maximum contrast, hence, in order to maximize the contrast range, it is possible to choose the sliding window size for a given pixel and speckle area and vice versa, it is possible to choose the ratio $M$ for a given sliding window size.

3.

Experimental Validation of Proposed Model

Laser light from a coherent source (Verdi, Coherent Inc., $\lambda =532\mathrm{nm}$) was expanded and collimated. Then the light passed through a diffuser to achieve a homogeneous illumination of a flow phantom, which consisted of a microfluidic slide (thinXXS Microtechnology AG) with channel diameter of $300\mu \mathrm{m}$. Intralipid (1%) was delivered with a syringe-based infusion pump into the channels at flow speeds from 4 to $20\mathrm{mm}/\mathrm{s}$ in step of $2\mathrm{mm}/\mathrm{s}$. Raw speckle were recorded with a CCD camera (Retiga 2000R) ($7.4\mu \mathrm{m}\times 7.4\mu \mathrm{m}$ pixel area). The camera was equipped with a macrolens with a variable aperture (Fig. 4) (see Ref. 38 for a detailed explanation of the experimental setup).

Fig. 4

Experimental setup used to test the spatial speckle correlation model. An expanded and collimated 532 nm laser illuminated the skin phantom. A Retiga CCD camera equipped with a variable-aperture macroacquired speckle images. A linear polarizer was placed in front of the lens to mitigate specular reflectance contributions to the images.

We modified the speckle size by changing the $f\#$ of the lens attached to the CCD camera:

Figure 5 shows the mathematical prediction of speckle contrast using Eq. (6) and the experimental data. The curves generated with Eq. (7) overlap when $M\ge 1$ but diverge when $M<1$. Considering no correlation between neighboring pixels, the speckle contrast $K_0$ calculated by the model approaches 1 when $M$ approaches zero (black continuous line) which is a contradiction, as the standard deviation should approach zero in that scenario and therefore the speckle contrast approaches zero, as observed experimentally. The speckle contrast must be $<1$ when $M$ is $<1$. When correlation subregions (i.e., $3\times 3$ and $5\times 5$) are considered, speckle contrast reduces for lower values of $M$ as expected.

Fig. 5

Solid lines show the contrast $K_S(3^2,1)$, $K_S(7^2,3)$, and $K_S({15}^2,)$ calculated using Eqs. (6) and (10) for spatial correlation submatrix of $3\times 3$, $7\times 7$, and $15\times 15\text{pixels}$ in pink, blue, and red, respectively. Dots show the experimental data for $3\times 3$, $7\times 7$, and $15\times 15\text{pixels}$, in pink, blue, and red, respectively.

The relationship between the pixel and the speckle grain size has an important effect on the speckle contrast.³⁹ Figure 5 shows spatial contrast calculated using a standard spatial contrast algorithm ($K_s(\mathrm{1,0})$), which does not consider the potential spatial correlation of adjacent pixels, and calculations with our analytic model considering a sliding window of $3\times 3$, $7\times 7$, and $15\times 15\text{pixels}$. For different window sizes, the curves have a different maximum contrast at different values of $M$. Therefore, with our model, knowledge of $M$ is sufficient to choose the appropriate sliding window and correlation subregion sizes to achieve the maximum speckle contrast.

We note that the experimental data in Fig. 5 agrees with calculations using our analytic model for all interrogated values of $M$. For $M\ll 1$ (the speckle grain size is larger than the pixel), choosing a submatrix at least as large as the speckle size to minimize loss of potential correlation between pixels is recommended.

Given a particular experimental setup and a sliding window size, $M$, $N$, and $p$ are defined. Therefore, the correction factor $\beta =K_s(N,p)$ [Eq. (6)] should be utilized to determine the perfusion index ($\mathrm{PI}=1/{\tau }_c$) from the contrast ($K$) obtained with the sliding window, and given by Eq. (2) or Eq. (4). Considering $T\gg {\tau }_c$, the PI is given by:³⁴ $\mathrm{PI}=\beta /T{K}^2$.

In Fig. 5, we use these window sizes merely as an example. However, the size of the sliding window is at the user’s discretion based on the characteristics of the speckle image they are analyzing. It should be noted that the larger the sliding window size is, the greater the loss of spatial resolution is.

Figure 6 illustrates the impact of the correction factor $\beta =K_s(N,p)$ on the PI image (for $T\gg {\tau }_c$), for two sliding window sizes (first row: $3\times 3$ and second row: $7\times 7\text{pixels}$) and a fixed value of $M=0.0338$. The first column, (a) and (d), shows the $\mathrm{PI}$ without considering the correlation between pixels [Eq. (5)], resulting in $\beta =K_s(\mathrm{1,0})=0.9827$. The second column, (b) and (e), shows the corrected $\mathrm{PI}$ [Eq. (6)] by considering all the potential correlations within the sliding window. In this case, the correction factors are $K_s(3^2,1)=0.7073$ and $K_s(7^2,3)=0.8198$, respectively. The third column, (c) and (f), is the difference between the first two columns. As shown, when $N=3^2$, the error in $\mathrm{PI}$ calculation is nearly 30% and when $N=7^2$, the error is $\sim 20\%$ with $M=0.0338$. The PI images shown in Fig. 6 represent the average PI over 30 frames.

Fig. 6

Accounting for the spatial correlation among pixels is required to calculate an accurate PI. (a)–(c) The $\mathrm{PI}=\beta /T{K}^2$ calculated using a sliding window of $3\times 3$. (d)–(f) A sliding window of $7\times 7$ pixel, with $M=0.0338$. (a), (d) The PI without considering the correlation between pixels, i.e., $\beta =K_s(\mathrm{1,0})=0.9827$. (b), (e) The corrected PI, $\beta =K_s(3^2,1)=0.7073$ and $\beta =K_s(7^2,3)=0.8198$, respectively. (c), (f) The difference between the first two columns.

From Fig. 6, we can observe that the difference between the first two columns is only the scalar multiplication by the correction factor, thus preventing the overestimation of the PI.

4.

Conclusions

We present a more accurate estimation of the constant $\beta$ and therefore of the contrast speckle contrast value by considering the effects of spatial correlation among the neighboring pixels. For a wide range of ratios of speckle grain size to pixel size, our model estimates the theoretical value of contrast more accurately than the standard model that does not consider the potential spatial correlation among adjacent pixels.³⁵ For the special case of Gaussian illumination, we present an analytical solution for the spatial speckle contrast of any spatial correlation submatrix size. The resultant impact of our new model is a more accurate estimation of speckle contrast, which in turn should lead to improved estimations of tissue blood flow.

5.

Appendix

The integrated intensity $I_{\eta ,\xi }$ over the physical area $a^2$ of the pixel $(\eta ,\xi )$ is given by

The mean intensity over the sliding window is

The expected value of the variance is

Given two pixels within the sliding window with coordinates $({\eta }^{\prime },{\xi }^{\prime })$ and $({\eta }^{″},{\xi }^{″})$, respectively, the parameter $1/{\mu }_{{\eta }^{″}-\eta \prime ,{\xi }^{″}-{\xi }^{\prime }}$ depends on the spatial correlation between the given pixels and $M$ the ratio between the pixel and speckle spot size. Due to it is not important, the actual pixel coordinates $({\eta }^{\prime },{\xi }^{\prime })$ and $({\eta }^{″},{\xi }^{″})$ but the distance between them $\eta ={\eta }^{″}-{\eta }^{\prime }$ and $\xi ={\xi }^{″}-{\xi }^{\prime }$. Therefore, the correlation factor $1/{\mu }_{\eta ,\xi }$ of the pixels $(\mathrm{0,0})$ and $(\eta ,\xi )$ it is identical to $1/{\mu }_{{\eta }^{″}-{\eta }^{\prime },{\xi }^{″}-{\xi }^{\prime }}$. The correlation factor is given by the following expression:

Computing $⟨i^2⟩$ from Eq. (14), where $⟨⟩$ denotes the ensemble average, and the summation spans all pixels within the sliding window of $\sqrt{N}\times \sqrt{N}$ pixels. Using the following notation:

It is possible to separate the summation into two sets of pixels (1) CW, where $(\eta ,\xi )\in (2p+1)\times (2p+1)$ and (2) OW, where $(\eta ,\xi )\notin (2p+1)\times (2p+1)$.

For set CW, we focus on particular cases. When $\eta =\xi =0$, we call this a central correlation (see Fig. 3). When $\eta =\xi \ne 0$, we call this a diagonal correlation. When $\eta \ne \xi =0$ or $\xi \ne \eta =0$, we call this a lateral correlation. Finally, when $(\eta ,\xi )\in (2p+1)\times (2p+1)$ and are not central, diagonal, or lateral correlations, we call these knight correlations due to the resemblance of the shape to the knight’s movement in chess. The other correlation, when the index $(\eta ,\xi )\in \mathrm{OW}$, is referred to as outsiders. Then Eq. (19) may rewritten as the sum of each type of correlations, i.e.:

For $p=0$, only the central correlation exists; hence, only the first summation in Eq. (20) is nonzero. Due to the symmetry in the correlation, some factors $1/{\mu }_{\eta ,\xi }$ are equivalent. For example, $1/{\mu }_{\eta ,\xi }$ is equivalent to $1/{\mu }_{\xi ,\eta }$, as well as any paired combination of $(\pm \eta ,\pm \xi )$. Approximately only the half of the pixels in one quadrant of the correlation subwindow are not equivalent to each other, without loss of generality, the fist quadrant is considered to classified/study these correlations, as it may see in Fig. 1(b). In this work, it has been classified the four different types of correlations: central, lateral, diagonal and knight.

For this reason, it is unnecessary to compute the correlation of all pixels. Instead, it is sufficient to compute the correlation of the pixels that are not equivalent and to weight their contributions appropriately.

Identifying the existence of $N$ central correlations, one for each pixel in the correlation subregion, is straightforward. The other types of correlation are not so straightforward, but it can be demonstrated that for the diagonal correlation type, $\eta =\xi \ge 1$, there are $4{(\sqrt{N}-\eta )}^2$ equivalent correlations. For lateral type, $\eta >\xi =0$, there are $4\sqrt{N}(\sqrt{N}-\eta )$ equivalent correlations. For the knight type, $1\le \xi \ne \eta \ge 2$, there are $8(\sqrt{N}-\eta )(\sqrt{N}-\xi )$ equivalent correlations and let $B(p)$ the number of outsider correlations.

For $p=0$, the correlation subwindow is only $1\times 1\text{pixel}$, i.e., only the autocorrelation is considered (central correlations), and Eq. (20) reduces to $⟨i^2⟩=\frac{1}{N^2}N⟨I_{\mathrm{0,0}}{I}_{\mathrm{0,0}}⟩$. For $p=1$, the correlation sub-window is $3\times 3\text{pixels}$, with only one central, one diagonal, and one lateral correlation. Therefore, Eq. (20) reduces to $⟨i^2⟩=\frac{1}{N^2}(N⟨I_{\mathrm{0,0}}{I}_{\mathrm{0,0}}⟩+4{(\sqrt{N}-1)}^2⟨I_{\mathrm{0,0}}{I}_{\mathrm{1,1}}⟩+4\sqrt{N}(\sqrt{N}-1)⟨I_{\mathrm{0,0}}{I}_{\mathrm{1,0}}⟩)$, which is the exact case analyzed by Skipetrov et al.³³ The $1\times 1$ pixels correlations subwindow ($p=0$) only has central correlation type, $3\times 3$ is the first subwindow ($p=1$) with lateral and diagonal correlations type, and the $5\times 5$ pixels correlation subwindow ($p=2$) is the first one that has one knight-type correlation, Eq. (20) reduce to

Coyotl-Ocelotl analyzed this case.³⁷

Henceforth $p\ge 2$, the first subregion with all types of correlation (diagonal, lateral, and knight) to avoid potential problems with missing terms in the summation. Then Eq. (20) is rewritten to

Another possible correlation is what we call outsider correlations, which expand beyond the total number of central, diagonal, lateral, and knight correlations within the correlation subregion $(2p+1)\times (2p+1)$:

Solving the summations, we obtain

Simplifying Eq. (24), we obtain

Here we assume that the outsider correlations are negligible by assuming that spatial correlations beyond the correlation subregion size are zero. For those pixels, $⟨I_{\mathrm{0,0}},I_{\eta ,\xi }⟩=⟨I_{\mathrm{0,0}}⟩⟨I_{\eta ,\xi }⟩={⟨I⟩}^2$, and using Eq. (18), we obtain

Using Eq. (23), some summation terms are eliminated with $B(p)$, and Eq. (27) becomes

Similarly, to compute $⟨c⟩$ [Eq. (16)], the term $⟨I_{\eta ,\xi }^2⟩=⟨I_{\eta ,\xi }{I}_{\eta ,\xi }⟩$ is equivalent to $⟨I_{\eta ,\xi }{I}_{\eta ,\xi }⟩=⟨I_{\mathrm{0,0}}{I}_{\mathrm{0,0}}⟩$. Using Eq. (18), Eq. (16) is rewritten as

After simplification:

Finally, dividing $⟨c⟩$ by ${⟨I⟩}^2$ results in the calculation of the spatial contrast: