Improved spatial speckle contrast model for tissue blood flow imaging: effects of spatial correlation among neighboring camera pixels

Abstract. Significance Speckle contrast analysis is the basis of laser speckle imaging (LSI), a simple, inexpensive, noninvasive technique used in various fields of medicine and engineering. A common application of LSI is the measurement of tissue blood flow. Accurate measurement of speckle contrast is essential to correctly measure blood flow. Variables, such as speckle grain size and camera pixel size, affect the speckle pattern and thus the speckle contrast. Aim We studied the effects of spatial correlation among adjacent camera pixels on the resulting speckle contrast values. Approach We derived a model that accounts for the potential correlation of intensity values in the common experimental situation where the speckle grain size is larger than the camera pixel size. In vitro phantom experiments were performed to test the model. Results Our spatial correlation model predicts that speckle contrast first increases, then decreases as the speckle grain size increases relative to the pixel size. This decreasing trend opposes what is observed with a standard speckle contrast model that does not consider spatial correlation. Experimental data are in good agreement with the predictions of our spatial correlation model. Conclusions We present a spatial correlation model that provides a more accurate measurement of speckle contrast, which should lead to improved accuracy in tissue blood flow measurements. The associated correlation factors only need to be calculated once, and open-source software is provided to assist with the calculation.


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][4][5] astronomy, [6][7][8] and ultrasound. 9,10iological 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,20Common uses of LSCI are to measure blood flow 21,22 and perform intraoperative assessment of tissue perfusion. 23or 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 24 and derived the following expression for the contrast K: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 1 ; 1 1 4 ; 5 2 0 where T is the camera exposure time, τ c is the correlation time of the backscattered light from the sample, σ is the standard deviation and hIi the average intensity.Since this seminal paper, the LSCI technique has evolved and therefore the contrast calculation too.For example, Lemieux and Durian 25 added a correction factor β (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, 26 resulting in g 2 ðτÞ ¼ 1 þ βjg 1 ðτÞj 2 , where g 1 ðτÞ is the electric field autocorrelation function and g 2 ðτÞ is the intensity autocorrelation function.Bandyopadhyay et al. 27 corrected an over simplification on the derivation of Eq. ( 1), applying a double integral of τ ¼ t 1 − t 2 , obtaining E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 2 ; 1 1 4 ; 3 7 8 where x ¼ T∕τ c .This expression still considers only the light scattered from dynamic scatterers and dimensionless pixels.Several groups 19,25,27,28 have since demonstrated the necessity of considering speckle not only from dynamic scatterers but also stationary scatterers.Static scatterers affect the Siegert relation 25 and thus the contrast, resulting in a revised form of g 2 : E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 3 ; 1 1 4 ; 2 9 4 where A ¼ and I s and I f represent the contributions from the static and dynamical scattered light, respectively: where ρ ¼ is the fraction of the fluctuating component of light that interacts with dynamic scatterers, and C n is an offset term that considers the noise contribution to the measurement.
Bandyopadhyay et al. 27 proposed a technique to measure β experimentally.Accounting for the finite size of the camera pixels 29 led to an analytical expression for β, assuming Gaussian illumination: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 5 ; 1 1 4 ; 1 3 4 where M is the ratio between the pixel area and the speckle area (i.e., the speckle grain size).β 1∕2 ðMÞ effectively corresponds to the spatial contrast with a sliding window of one pixel size.Equation ( 5) describes the effect of pixel and speckle sizes on the theoretical calculation of speckle contrast.However, it does not consider the likely correlation among neighboring pixels.
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 ≥ 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Þ.

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 β and therefore of the contrast.
Skipetrov et al. 33 employed a sliding window of 3 × 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ðrÞ recorded by a digital camera, where r ¼ x î þ y ĵ is the position vector on the CCD image.Here we consider a sliding window of ffiffiffiffi N p × ffiffiffiffi N p pixels, with N being an odd number.We define the speckle spot size as comprising a subregion of . 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). 2 g 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; ηÞ, ðη; 0Þ, ð0; −ηÞ, and ð−η; 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 ffiffiffiffi N p × ffiffiffiffi N p pixels, and hence an arbitrary correlation subregion of ð2p þ 1Þ × ð2p þ 1Þ pixels may be written as (details can be found in Appendix): E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 0 6 ; 1 1 7 ; 2 3 3 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 ≥ 2, i.e., a 5 × 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 ≥ 2. 1∕μ η;ξ is the correlation factor of the pixel with coordinates ðη; ξÞ and with the central pixel ð0;0Þ of the sliding window, given by where g 1 is a first-order correlation function.At this point, we have not considered the correlation function ðg 2 1 ðr − r 0 ÞÞ; this will be the topic of the next section.The factor β in Eq. ( 4) corresponds to the spatial contrast K 2 s ðN; pÞ (i.e., β ¼ K 2 s ðN; pÞ) for any window size, while Eq. ( 5), K 2 s ð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.
As it seen in Table 1, the contribution of each lateral, diagonal, and knight correlation decreases with increasing subregion size.

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

Analytical solutions
The Gaussian correlation function g i 33,34 may be written using ; t e m p : i n t r a l i n k -; e 0 0 8 ; 1 1 7 ; 6 4 4 where A c is the correlation area (i.e., speckle size).For a speckle area with radius b: ; t e m p : i n t r a l i n k -; e 0 0 9 ; 1 1 7 ; 5 9 4 Substituting Eq. ( 9) into Eq.( 7), we obtain ; t e m p : i n t r a l i n k -; e 0 1 0 ; 1 1 7 ; 5 4 3 where M ¼ a 2 ∕π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: 29 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 1 ; 1 1 7 ; 3 6 8 Equation ( 11) is identical to the correction factor β [see Eq. ( 5)] proposed in Refs.19 and 25.

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).
The parameters 1∕μ 0;0 , 1∕μ 1;0 , and 1∕μ 1;1 are sufficient to calculate the contrast K for a spatial correlation subregion of 3 × 3 pixels. 36For a subregion of 5 × 5 pixels, we also should include the terms 1∕μ 2;0 , 1∕μ 2;1 , and 1∕μ 2;2 . 37Now it is possible to calculate the speckle contrast using Eq. ( 10) for any arbitrary odd size of sliding window ( ffiffiffiffi N p ) and subregion (2p þ 1). Figure 3 shows the plots of the analytical (solid lines) and numerical (symbols) contrast for a correlation subregion of 1 × 1 (i.e., no spatial correlation), 3 × 3, and 5 × 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 σ 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.10) needed to calculate the speckle contrast of a 7 × 7 sliding window.Central Fig. 2 Comparison between the analytical expression (solid lines) using Eq. ( 10) and a numerical solution (symbols) of Eq. ( 7) using Mathematica 11.
In addition, the integration of multiple grains per pixel leads to a reduction in intensity variance in the speckle pattern.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.

Experimental Validation of Proposed Model
Laser light from a coherent source (Verdi, Coherent Inc., λ ¼ 532 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 μm.Intralipid (1%) was delivered with a syringe-based infusion pump into the channels at flow speeds from 4 to 20 mm∕s in step of 2 mm∕s.Raw speckle were recorded with a CCD camera (Retiga 2000R) (7.4 μm × 7.4 μ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).
We modified the speckle size by changing the f# of the lens attached to the CCD camera: where Mag is the magnification of the lens.The spatial speckle contrast was calculated from the experimental data using a sliding window of N ¼ 3 2 ; 7 2 , and 15 2 .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 ≥ 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 × 3 and 5 × 5) are considered, speckle contrast reduces for lower values of M as expected.
The relationship between the pixel and the speckle grain size has an important effect on the speckle contrast. 39Figure 5 shows spatial contrast calculated using a standard spatial contrast algorithm (K s ð1;0Þ), which does not consider the potential spatial correlation of adjacent pixels, Fig. 3 Solid lines show the contrast calculated using Eq. ( 6) for a correlation subregion of 1 × 1 (no spatial correlation), 3 × 3, and 5 × 5 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 × 1, 3 × 3, and 5 × 5 pixels.Good agreement is observed between the numerical and analytical expressions.
and calculations with our analytic model considering a sliding window of 3 × 3, 7 × 7, and 15 × 15 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 ≪ 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 β ¼ K s ðN; pÞ [Eq.( 6)] should be utilized to determine the perfusion index (PI ¼ 1∕τ c ) from the contrast (K) obtained with the sliding window, and given by Eq. ( 2) or Eq. ( 4).Considering T ≫ τ c , the PI is given by: 34 PI ¼ β∕TK 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.
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 × 3, 7 × 7, and 15 × 15 pixels in pink, blue, and red, respectively.Dots show the experimental data for 3 × 3, 7 × 7, and 15 × 15 pixels, in pink, blue, and red, respectively.(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 PI calculation is nearly 30% and when N ¼ 7 2 , the error is ∼20% with M ¼ 0.0338.The PI images shown in Fig. 6 represent the average PI over 30 frames.
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.

Conclusions
We present a more accurate estimation of the constant β 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. 35For 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.

Appendix
The integrated intensity I η;ξ over the physical area a 2 of the pixel ðη; ξÞ is given by E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 3 ; 1 1 7 ; 1 1 4 The mean intensity over the sliding window is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 4 ; 1 1 4 ; 7 2 4 and the corresponding unbiased variance is E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 5 ; 1 1 4 ; 6 7 2 The expected value of the variance is ; t e m p : i n t r a l i n k -; e 0 1 6 ; 1 1 4 ; 6 1 9 where h i denotes the ensemble average.Given two pixels within the sliding window with coordinates ðη 0 ; ξ 0 Þ and ðη 00 ; ξ 00 Þ, respectively, the parameter 1∕μ η 00 −η 0;ξ 00 −ξ 0 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 ðη 0 ; ξ 0 Þ and ðη 00 ; ξ 00 Þ but the distance between them η ¼ η 00 − η 0 and ξ ¼ ξ 00 − ξ 0 .Therefore, the correlation factor 1∕μ η;ξ of the pixels ð0;0Þ and ðη; ξÞ it is identical to 1∕μ η 00 −η 0 ;ξ 00 −ξ 0 .The correlation factor is given by the following expression: using Eqs.( 13) and ( 17), we can rewrite hI 0;0 I η;ξ i as ; t e m p : i n t r a l i n k -; e 0 1 8 ; 1 1 4 ; 4 4 2 Computing hi 2 i from Eq. ( 14), where h i denotes the ensemble average, and the summation spans all pixels within the sliding window of ffiffiffiffi N p × ffiffiffiffi N p pixels.Using the following notation: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 1 9 ; 1 1 4 ; 3 7 7 It is possible to separate the summation into two sets of pixels (1) CW, where ðη; ξÞ ∈ ð2p þ 1Þ × ð2p þ 1Þ and (2) OW, where ðη; ξÞ ∈ = ð2p þ 1Þ × ð2p þ 1Þ.
For set CW, we focus on particular cases.When η ¼ ξ ¼ 0, we call this a central correlation (see Fig. 3).When η ¼ ξ ≠ 0, we call this a diagonal correlation.When η ≠ ξ ¼ 0 or ξ ≠ η ¼ 0, we call this a lateral correlation.Finally, when ðη; ξÞ ∈ ð2p þ 1Þ × ð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 ðη; ξÞ ∈ OW, is referred to as outsiders.Then Eq. ( 19) may rewritten as the sum of each type of correlations, i.e.: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 0 ; 1 1 4 ; 2 3 3 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∕μ η;ξ are equivalent.For example, 1∕μ η;ξ is equivalent to 1∕μ ξ;η , as well as any paired combination of ðAEη; AEξÞ.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, η ¼ ξ ≥ 1, there are 4ð For p ¼ 0, the correlation subwindow is only 1 × 1 pixel, i.e., only the autocorrelation is considered (central correlations), and Eq. ( 20) reduces to hi 2 i ¼ 1 N 2 NhI 0;0 I 0;0 i.For p ¼ 1, the correlation sub-window is 3 × 3 pixels, with only one central, one diagonal, and one lateral correlation.Therefore, Eq. ( 20) reduces to hi 0;0 I 1;0 iÞ, which is the exact case analyzed by Skipetrov et al. 33 The 1 × 1 pixels correlations subwindow (p ¼ 0) only has central correlation type, 3 × 3 is the first subwindow (p ¼ 1) with lateral and diagonal correlations type, and the 5 × 5 pixels correlation subwindow (p ¼ 2) is the first one that has one knight-type correlation, Eq. ( 20) reduce to E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 1 ; 1 1 7 ; 5 2 1 0;0 I 2;1 iÞ: Coyotl-Ocelotl analyzed this case. 37enceforth p ≥ 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 E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 2 ; 1 1 7 ; 4 0 3 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Þ × ð2p þ 1Þ: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 3 ; 1 1 7 ; 2 8 7 Solving the summations, we obtain E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 4 ; 1 1 7 ; 2 1 7 Simplifying Eq. ( 24), we obtain E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 2 5 ; 1 1 7 ; 1 4 2 Here we assume that the outsider correlations are negligible by assuming that spatial correlations beyond the correlation subregion size are zero.For those pixels, hI 0;0 ; I η;ξ i ¼ hI 0;0 ihI η;ξ i ¼ hIi 2 , and using Eq. ( 18), we obtain Equation ( 33) is an analytic expression for speckle contrast that assumes both an arbitrary sliding window size and an arbitrary correlation subregion of pixels (p ≥ 2).At this time, we do not consider any heterogeneity in illumination intensity ðg 2 1 ðr − r 0 ÞÞ; this should be the subject of future work.
Equation ( 33) is expressed in terms of the factors 1∕μ η;ξ .Those factors are given by Eq. ( 17), which may be simplified using the following changes of variables w ¼ x − x 0 , w ¼ ðx þ x 0 Þ∕2, u ¼ y − y 0 , and u ¼ ðy þ y 0 Þ∕2 equivalent to the triangular weight reported in the literature after Bandyopadhyay et al. 27 Equation ( 17) can thus be rewritten as To this point, we have not made any assumptions about the shape of the correlation function g 2 1 ðw; uÞ.We now assume a Gaussian correlation shape: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 5 ; 1 1 7 ; 3 9 4 where A c is the correlation area and has been approximated to a circle of radius b.Substituting Eq. ( 35) into Eq.(34), we arrive at the following expression for the factors 1∕μ η;ξ , which is Eq.(10) in the main text: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 3 6 ; 1 1 7 ; 3 2 0

Disclosures
No conflicts of interest, financial or otherwise, are declared by the authors.

Fig. 1
Fig. 1 (a) Schematic representation of the sliding window of ( ffiffiffiffi N p × ffiffiffiffi N p ) pixels (pink) and the correlation subregion of ð2p þ 1Þ 2 (blue).In this example, the sliding window of size ffiffiffiffi N p ¼ 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).

Fig. 4
Fig.4Experimental 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.

Figure 6
Figure 6 illustrates the impact of the correction factor β ¼ K s ðN; pÞ on the PI image (for T ≫ τ c ), for two sliding window sizes (first row: 3 × 3 and second row: 7 × 7 pixels) and a fixed value of M ¼ 0.0338.The first column, (a) and (d), shows the PI without considering the correlation between pixels [Eq.(5)], resulting in β ¼ K s ð1;0Þ ¼ 0.9827.The second column, (b) and (e), shows the corrected 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 PI calculation is nearly 30% and when N ¼ 7 2 , the error is ∼20% with M ¼ 0.0338.The PI images shown in Fig.6represent the average PI over 30 frames.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.

Finally, dividing hci by hIi 2
results in the calculation of the spatial contrast: E Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0

E
Q -T A R G E T ; t e m p : i n t r a l i n k -; e 0 uÞðaðη − 1Þ þ wÞdw du:

Table 2
Correlation factors calculated using Eq. (