Laser speckle contrast imaging (LSCI) is a promising technique for the non-invasive measurement of dynamic systems (i.e., flowmetry), such as blood flow [1–6]. When coherent light scatters off of particles, a speckle image is formed. When those particles are in motion, the speckles will be dynamic as well. The result is both speckle translation and speckle boiling [7, 8], with a characteristic speckle decorrelation time τ_c. When imaged with a camera with exposure time T, the resulting speckle image will undergo blurring, the amount of which depends on the particles’ velocities. The blurring is quantified using the speckle contrast, K, thus providing a metric for velocity, V [9]. At present, this makes LSCI useful for relative velocity mapping [1]. Due to its simplicity, high spatial resolution, and low cost, LSCI has already been widely adopted [10–14].

However, quantitative measurements with LSCI have been elusive due to the lack of a method for quantitatively determining the velocity from K, which prevents us from making quantitative measurements with LSCI. Another optical technique for velocimetry is particle image velocimetry (PIV), which, although established, has the disadvantage that direct imaging is required and is thus invasive in nature. In a recent paper, the velocity profile was quantitatively reconstructed with the new optical speckle image velocimetry technique [15] that combines LSCI with PIV, but their technique relies on negligible speckle boiling, which is only attainable with invasive measurements. The sidestream dark field LSCI technique [16] still requires some direct imaging of the flow as well. Efforts have been made to overcome these problems for non-invasive measurements, such as improvements to the analytical relationship between K and τ_c [17, 18], studying the effect of the temporal correlation function of light [5, 19, 20], and using multiple exposure times [21, 22] or multiple wavelengths [23]. Although LSCI is promising, much work is still required to make it a fully quantitative measurement technique.

In this work, we continue the work of Duncan and Kirkpatrick [17] by investigating how K depends on τ_c and τ_c in turn on V. Once we accurately know these relationships and what they depend on, we can convert a measured K into the velocity of the scattering system, thus enabling us to make quantitative measurements with LSCI. To obtain those relationships, we perform computer simulations using our new in-house code [24].

The code is based on Mie theory, which describes the scattering of a linearly polarised plane wave by a single homogeneous spherical particle. Using a far-field approximation, the plane wave that was scattered by a particle locally becomes a plane wave again. In that manner, multiple scattering between particles is implemented iteratively, in which each particle scatters to each other particle, including backscattering, until successive scattering orders become negligible. Finally, all scattered fields are gathered at a two-dimensional grid of infinitesimal points (i.e., our “simulated camera”), at which the intensity I is calculated. Using a separate computational fluid dynamics code, we then evolve the particles in time. The instantaneous light scattering calculation is repeated at n_int rapid time intervals and then averaged over to mimic the finite integration time T of a camera. The result is a fully interferometric code, capable of simulating dynamic speckle.

We do not simulate any kind of imaging system such as lenses. Consequently, we study objective speckle as opposed to the subjective speckle that forms in the imaging plane of a lens. Since both types of speckle have similar dynamics, simulating lenses is irrelevant to our simulation.

Several simplifications were made that enable us to study dynamic systems within a reasonable computational time. The strongest is the aforementioned far-field approximation in Mie theory. This is easily violated, as this requires an interparticle distance of δr ≫ 0.1 mm (i.e., for our parameters outlined below). The second part of the far-field approximation is that the particle size should be ≪ δr, which is much more easily satisfied than the previous assumption. However, although these assumptions are not valid for, e.g., real blood flow (δr ~ 10⁻⁵ m), they are still satisfied in our simulations, because we are limited to relatively few particles by computational constraints. Thus strictly speaking our model is only applicable to sufficiently dilute flow, but it may be expected that the results on dynamic speckle imaging are more widely applicable nonetheless. More details about our code and the assumptions may be found in our previous paper [24].

To study how K ≡ σ_I/〈I〉 depends on V, we use a simple cylindrical geometry with plug flow (i.e., a uniform constant velocity profile). The cylinder is 1 cm long with a 1 mm radius, which is characteristic for the external carotid artery.1 The camera and the laser are placed at right angles, and also orthogonal to the cylinder’s axis (see Fig. 1). Hundred tracer particles with radius 4 m are randomly distributed over the geometry. When a particle leaves the cylinder, it is reinserted at the entrance at the same radial and polar position (i.e., cyclic boundary conditions). Although our code is capable of simulating more complex geometries and flow profiles, we chose this simple setup as a first step to minimize the effect of speckle boiling. Effects of using a different flow profile on the noise induced by the associated increase in speckle boiling have been studied in our previous paper [24].

On the optics side of the code, we use a real refractive index of 1.52 for the spheres and 1.00 for the surrounding medium. The wavelength of illumination is 532 nm. The camera is placed at a distance of 25 cm from the cylinder, and is given a size of 1.25 cm by 5.00 cm with 128² pixels.2 With these settings speckles are underresolved, on average using merely 1 pixel for every 25 speckles (i.e., 0.2² pixels per typical speckle). Whereas an experiment should satisfy the Nyquist criterion as to prevent an artificial reduction of K due to spatial averaging over the finite size of a pixel [26], our simulation uses infinitesimal point pixels and thus does not have this problem. Consequently, we can use fewer pixels that are separated by multiple speckle sizes to obtain a better statistical representation of the speckle space while retaining all intensity fluctuations [24]. Finally, from the obtained speckle pattern, the speckle contrast is calculated using local speckle contrast analysis [20, 24] with windows of 8 × 8 pixels.

A study of the effect of velocity V and camera integration time T on the speckle contrast K is shown in Figure 2. The simulations were performed with 10 different sets of random initial particle positions, which allowed us to calculate the systematic error caused by fluctuations in the precise particle positions that will always be present in experiments and simulations alike.3 The data points in the figure are the resulting mean speckle contrasts, and the error bars show the resulting standard deviations. Additional errors are not incorporated into these error bars, such as due to spatial and temporal discretization (which are ~2% each [24]). The figure nicely shows that K ≈ 1 at zero velocity, which is the expected value for fully-developed static speckle [27]. Actually, K is slightly (~1.2%) less than unity [24], but this is not substantial when compared to the uncertainties of discretization and of the random particle instantiations.

As either V or T is increased, K decreases because more motion is captured. Generally speaking, the error is also lower for higher V and T, presumably because for those situations averaging occurs over a larger range of particle positions, which renders the random particle instantiations less relevant.

We hypothesize that the relevant parameter should in fact be d = VT, which is the physical distance the particles have travelled during the camera integration time. Figure 3 shows the same data as Figure 2, but rescaled using d. Clearly, all data points nicely collapse onto a single master curve. In fact, in a simulation with the same initial particle instantiations, doubling either V or T gives precisely identical particle positions and thus identical speckle and K.

Several analytical expressions for K are already well-known [17, 19, 20]. As a numerical experiment, we compare our simulations with these analytical models as to confirm/determine their applicability. Whereas physical experiments are good at measuring at large integration times, our simulations can be used for very small integration times not yet approachable in experiments. More generally, simulations can be used for circumstances that are difficult to reach experimentally, and thereby compliment experimental results.

The analytical expressions are derived from the temporal fluctuation statistics of the speckles caused by the motion of the scatterers. Using the autocovariance of the temporal fluctuations $C_t^{\left(2\right)}\left(\tau \right)$, it is possible to integrate over a time period T to obtain the spatial variance σ_I of time-integrated speckle [20, 28], and thus also K:(1)$${\sigma }_I^2\left(T\right)=\frac{2}{T}{\int }_0^T\mathrm{}\left(1-\frac{\tau }{T}\right)C_t^{\left(2\right)}\left(\tau \right)\mathrm{d}\tau .$$

However, the difficulty in deriving an analytical expression is that the motion of the scatterers is usually a priori unknown or too complicated; therefore, $C_t^{\left(2\right)}\left(\tau \right)$ caused by said motion must be assumed.

When a Lorentzian covariance is assumed, the speckle contrast is:(2)$$K^2=\beta \frac{{\tau }_c}{2T}\left[2-\frac{{\tau }_c}{T}\left(1-e^{-2T/{\tau }_c}\right)\right],$$ where β is a factor that corrects for the loss of correlation due to the ratio of pixel size to speckle size [20], with β = 1 for infinitesimal pixels. Analogously, for an assumed Gaussian covariance, it holds that:(3)$$K^2=\beta \frac{{\tau }_c}{2T}\left[\sqrt{2\pi } \mathrm{erf}\left(\frac{\sqrt{2}T}{{\tau }_c}\right)-\frac{{\tau }_c}{T}\left(1-e^{-2(T/{\tau }_c{)}^2}\right)\right].$$

The Lorentzian equation is valid for unordered (Brownian) motion, whereas the Gaussian equation is more appropriate for ordered motion. Since those two types of motion are statistically independent, reality is likely somewhere between these two limits [19].

However, it is not yet clear whether these relationships are truly applicable in practice. An arguably more appropriate relationship for blood flow was derived while assuming a constant velocity [29]:(4)$$K^2=\beta \frac{{\tau }_c}{3.83T}{\int }_0^{3.83T/{\tau }_c}\mathrm{}{\left(\frac{2J_1\left(x\right)}{x}\right)}^2\mathrm{d}x,$$ where J₁(x) is the Bessel function of the first kind. The number 3.83 was introduced to define τ_c in equation (1) as the time τ after which J₁ (x) = 0 for the first time. This curve follows the Gaussian curve closely for small T/τ_c, whereas it follows the Lorentzian for large T/τ_c. This is sensible, because for small T/τ_c the motion is very ordered (~Gaussian), whereas for large T/τ_c the blurring caused by the accumulated speckle blurring is so large that we might as well have been looking at speckle boiling caused by unordered motion (~Lorentzian).

These expressions all give K as a function of τ_c, whereas the goal is to measure V. To that end, an expression for τ_c(V) is needed. It was postulated that τ_c scales inversely with V: [17](5)$${\tau }_c=\frac{w}{V},$$ where the proportionality constant w should be a characteristic length scale; however, its value is still being disagreed upon by two orders of magnitude [5]. The first proposal simply used the wavelength, w_λ = λ/2π, although without any derivation [1, 17]. Recent experimental research on speckle dynamics still uses this relation [30]. Later it was suggested that w should be the speckle size [17, 19],(6)$$w_{\mathrm{speckle}}=\frac{\lambda z}{D},$$ where z is the distance between the object and image planes, and D is the imaging aperture.4 This expression makes sense physically, because in the time period w_speckle/V the speckles have first translated a relevant distance to start decorrelating, which is precisely what the time τ_c describes. Since there are no lenses in our numerical setup, our “aperture” (D) is the illuminated area, which is our entire cylinder.

We may compare our results with these expressions by first noting that all expressions for K in equations (2)–(4) depend only on the ratio τ_c/T. Upon substituting equation (5), we see that the analytical results for K all depend on the quantity w/VT. This is consistent with our independent finding that d = VT is the relevant parameter for studying K(V), as was evidenced by Figure 3. We view this finding as evidence that τ_c does indeed scale inversely with V, cf. equation (5).

Having written K(τ_c/T) as K(w/d), equations (2)–(4) may be used to fit our simulation results with fitting parameters w and β. Figure 4 shows the same data as shown in Figures 2–3, but now with the x-axis rescaled to (τ_c/T)/w = d⁻¹. This scale was chosen to resemble the usual τ_c/T scale [17, 20], but deviates because we used w to fit our data with known V. Consequently, τ_c in equation (5) is not well-defined, as it differs for each fit (by up to a factor two). For the fit of equation (4) we have used the exact solution of the integral. For the fit of the Lorentzian model β = 1 is used, because the fit would result in β > 1, resulting in non-physical K > 1 at low velocities. Finally, all three models would have been bad if we had fit them across the whole domain; therefore, only the region of large τ_c/T contributes to the fit.

From the figure, we note that the Lorentzian model describes our results poorly, whereas the other models do an excellent job for large τ_c/T ~ d⁻¹ (low V). This is consistent with the fact that we have flow at a constant velocity, whereas the Lorentzian model is more appropriate for Brownian motion [17]. The Gaussian and constant-velocity models resemble each other more closely, as they are both for ordered motion.

However, β describes the reduction of K due to the finite pixel size, whereas we have infinitesimal pixels. Therefore, we may not use β as a fitting parameter, and actually should just take β = 1, as is shown in Figure 5. Note that β has a negligible influence at small τ_c/T ~ d⁻¹ (high V), and mostly results in a vertical shift at large d⁻¹ (low V). The main difference here is that the analytical models have their asymptote at K = 1, whereas our simulations yield a value slightly below one. The reason is that our simulations implement multiple scattering, which reduces K. However, in our rather dilute simulations, single scattering still contributes significantly (~90%), resulting in only a slightly lower value for K. The analytical models do not incorporate this effect, and thus β provides a first order approximation to including multiple scattering analytically.

Next, it may be seen that the Gaussian model performs slightly better across the whole domain than the constant-velocity model, although our simulations do have a constant velocity. A possible explanation could have been that particles move in and out of the laser’s view, which results in a small amount of speckle boiling. However, simulations without cyclic boundary conditions (in which we have a purely translational speckle) have revealed a less than 1% difference in the results; thus speckle boiling is not the cause.

Related to that, all models describe the results at small τ_c/T ~ d⁻¹ (high V) poorly, as K seems to saturate in the simulations, which is why the fits were made for large d⁻¹. However, fitting the models at small d⁻¹ instead does not yield a good fit (not shown); therefore, the models cannot describe the simulations in this regime. In an experiment, the effect of static scatterers would be to increase the minimum contrast value [20], which is precisely the effect we observe. However, the simulations do not have any static scatterers. Hypotheses for the difference are effects of single versus multiple scattering, and related thereto the diluteness of our flow. More specifically, the theoretical models do not incorporate multiple scattering, but its contribution to the simulated intensity is still minor in the present simulations regardless. Therefore, a more likely explanation is the diluteness of our simulations in combination with the major contribution from singly scattered light. The combination causes interferometric fringes to appear on our camera. Although local contrast analysis (i.e., windowing) diminishes their influence, they do have the effect of increasing K and are not affected by the blurring at high velocities (low d⁻¹) [24]. Alternatively, it could also well be that the theoretical models are inadequate in this regime, as they, too, make assumptions on the scattered field’s statistics. Future research will need to point out which of the above hypotheses is the case. Nevertheless, our results are meaningful at large d⁻¹, which is the regime that is most difficult to study experimentally due to the requirement of a camera with small integration time.

Finally, it is of prime interest to study the obtained value for the fitted w, as there exists disagreement in the literature. The hypothesised expressions for w in and just above of equation (6) yield: w_λ = 0.0847 m and w_speckle = 13.3 m × 65.5 m. w_speckle is given two values, as our aperture5 is rectangular with aspect ratio 5, and thus our speckles, too, are rectangular with aspect ratio 1/5 [24]. However, in the direction of motion the speckles have width w_speckle = 13.3 m, which is the only relevant length scale for decorrelating speckles due to translational speckle in that direction.

Table 1 shows the obtained values of w and β from the fits in Figures 4 and 5, with which we can compare w_λ and w_speckle. It is clear that w_λ is off by several orders of magnitude, and thus is inadequate. w_speckle, on the other hand, is strikingly close to our results, and in particular to the Gaussian, being also the best fit. Thus w ≈ 1.06 w_speckle in combination with the Gaussian model best describes our simulation results.6 Although this does not proof the correctness of equation (6), for which future research should investigate the influence of λ, z and D, it does make w_speckle an extremely likely candidate.

In summary, we have presented results of our new computer code, which simulates how a plane wave of coherent light scatters off of a collection of moving particles using Mie theory to form a dynamic interferometric speckle pattern (i.e., Laser Speckle Imaging (LSI)). By mimicking the finite integration time T of a real camera, we have shown how the speckle contrast K depends on particle velocity V and on T. Existing theoretical models already describe how K depends on the speckle decorrelation time τ_c and on T, and it is believed that τ_c = w/V; although the value of w is still disagreed upon in the literature. We provide evidence that τ_c does indeed scale inversely with V, and that w_speckle = λz/D (multiplied by an O(1) constant) is a very likely candidate for w. The Gaussian correlation model does an excellent job at describing our simulation results for large τ_c/T (low V), but deviates considerably for low τ_c/T (high V), for which we provide several hypotheses for future research. However, the Lorentzian model is unsuitable for ordered flow (i.e., advection). Other optical scattering techniques using photon correlation methodologies also have the continuing discussion about decomposing flow into advection and Brownian motion [31].

The strength of simulations, once validated, is that circumstances that are difficult to reach experimentally may be studied using numerical experiments. Therefore, in future research our computer code may be used to study the effect of all relevant parameters – and in particular those that are not easily accessible in experiments – which will help develop LSI as a fully quantitative non-invasive measurement technique for flowmetry in turbid media (e.g., blood flow7).