1. Introduction

An accurate and often continuous assessment of microvascular, regional blood flow has many implications for diagnosis and treatment of diseases and for the study of healthy physiology. Despite continued efforts to establish practical means for measuring microvascular, regional blood flow in a non-invasive manner, this remains an important unmet need. One potential approach uses near-infrared, coherent light and the arising speckles after its diffusion [1–4].

Coherent laser light can be used to non-invasively measure local microvascular blood flow in tissue by detecting the fluctuating speckle patterns after light interaction with the tissue [5–9]. For the purposes of this manuscript, we will focus on deep-tissue, i.e. those that utilize light that penetrates up to several centimeters, measurements using photon diffusion. This is possible since near-infrared ($\sim$650-1000 nm) light is only mildly absorbed in most tissues.

In the field of near-infrared diffuse optics, there are two common methods for determining blood flow from laser speckles. The first consists of measuring the speckle intensity autocorrelation ($g_2(\tau )$) or the electric field autocorrelation ($g_1(\tau )$) over a continuous range of decay times ($\tau$) to derive a blood flow index [10]. Diffuse correlation spectroscopy (DCS) [10–12] and its variants [13–15] utilize this method for quantifying the speckle statistics to determine blood flow. The second common method consists of quantifying the speckle intensity statistics using a parameter called the “speckle contrast” ($\kappa$). Several related techniques measure $\kappa$ to measure blood flow. These include tomographic techniques (SCOT, scDCT) for the three-dimensional mapping of blood flow from measurement of $\kappa$ [16,17] and techniques to calculate one or two-dimensional maps of blood flow (DSCA, SCOS, LSF, LASCA, LSCI, iSVS) [2,8,18–22]. Of these, some techniques (LASCA and LSCI) are non-diffuse methods and therefore only measure superficial blood flow [8,20].

Diffuse optical methods using the laser speckle contrast can achieve similar blood flow information as DCS at an overall cheaper cost per detector channel since $\kappa$ is an integral of $g_2(\tau )$ over the delay times up to a longer exposure time. In other words, common scientific cameras can be utilized as “slower” detectors. This claim has been supported by experiments [3,23–25], simulations [26], and most recently by theoretical analyses [27].

A thorough analysis of the measurements utilizing the intensity auto-correlation of the speckle statistics, i.e. DCS, has previously been developed and tested [28–33]. Among other uses, these works have allowed the design of components (detectors, sources) and systems that target specific goals in detection precision and accuracy in DCS.

Despite the increasing prevalence in literature of the use of speckle contrast techniques, a comprehensive method for determining the effects experimental parameters have on the accuracy and precision of $\kappa$ has not yet been developed. Accuracy in speckle contrast values, particularly in scenarios with low levels of detected light, is important to consider as the effects of detector noise can greatly influence the detected signal. Valdes et al. [2] first described this phenomenon, and subsequently developed a noise removal algorithm to reduce the effect of detector noise on the measured value of $\kappa$. This algorithm has been shown to be effective, however it does not correct for all detector effects, in particular shot noise.

Previous work to optimize accuracy and precision in speckle contrast measurements includes theoretical and experimental characterization of the sampling of speckles on the precision of measured $\kappa$ [34–36], and the effect of the imaged speckle to camera pixel ratio on the accuracy of $\kappa$ [37–39]. These earlier works did not account for the effect of experimental sources of noise, particularly detector noise, on the measured accuracy and precision of the speckle contrast signal. Recently, this gap in the existing literature was addressed by Zilpelwar et.al. [40] through a simulation method which modeled the generation and detection of decorrelating speckles including detector noise effects. The authors demonstrated that the developed model is able to simulate both the values of $\kappa$ as well as the noise in $\kappa$ detected using sCMOS cameras. Using this simulation, the authors investigate the effect of speckle to pixel size ratio, exposure time, and detected photon count rate on $\kappa$ and its signal to noise ratio (SNR) for two commercially available cameras.

We have developed a separate simulation model to Zilpelwar et.al. [40], but with a similar aim of simulating the behavior of $\kappa$ with respect to detector noise and other experimental parameters. Our model addresses details not included in Ref. [40] such as the efficacy of the detector noise correction by Valdes et.al. [2], and the behavior of $\kappa$ in a multi-scattering regime in a semi-infinite geometry. We are specifically interested in characterizing the accuracy and precision of speckle contrast measurements taking into consideration experimentally relevant parameters such as the noise specifications of the detectors, the exposure time of the experiments, the detected photon-count rate, the measured medium, and the sampling of the detected speckles. To this end, the developed method was first verified experimentally for its ability to simulate $\kappa$ and the noise in $\kappa$. After verifying the simulation method, the method was used to study the effect of accuracy and precision of $\kappa$ in various experimental scenarios. Finally, the simulations were used to design and optimize a system capable of measuring baseline cerebral blood flow non-invasively in an adult human.

2. Methods

Here we focus on two dimensional detectors ($i \times j$) with “pixels” but the results can be generalized to other standard detectors. As will be evident later on, it is more convenient to use the square of the speckle contrast ($\kappa ^2$) for the analysis. We assume that the $\kappa ^2$ is derived from sampling $n$ speckles that are distributed over space ($w_z$) and/or over time by repeated measurements ($w_t$). These $n$ speckles sampled over $w_z$ and/or $w_t$ are used to estimate the probability distribution of the speckle intensity. From these $n$ speckles, the mean intensity ($\mu (I)$) and the variance of intensity ($\sigma ^2(I)$) are determined.

Even in the case of ideal detectors and light sources, the calculated values are not exactly equal to the true mean and the true variance due to the effects of limited sampling. In experiments, the situation is more complex due to additional sources that contribute to the observed photon statistics such as the detector noise which further influence the measured values of mean and variance.

Therefore, these measurement effects must be accounted for in order to experimentally determine a “corrected $\kappa ^2$”, or the best estimate of the true value of $\kappa ^2$. For common detectors, these corrections include a dark frame subtraction which attempts to account for the dark and read-out signal and a statistical correction attempting to estimate the shot noise as well as the dark and read-out noise variances [2].

The speckle contrast is an alternative data-type that is used to characterize the decorrelation time ($\tau _c$) of the intensity autocorrelation of the speckle statistics which is more commonly utilized [27,41]. $\tau _c$ is in turn dependent on several aspects such as the the optical properties of the medium, the dynamics of the scatterers, the measurement geometry, the source wavelength and more. The signals that are detected in a common detector are affected by this statistical profile which in turn affects the noise statistics. Therefore, in order to simulate realistic speckle contrast signals, we need to take all this into account and incorporate the appropriate aspects of the detectors. An illustrative flowchart of the method that has been developed is shown in Fig. 1 and is further detailed below.

 

Fig. 1. Flow chart for simulating frames of correlated speckles and $\kappa ^2$. These simulations aim to simulate a variety of experimental setups such as in sub-figure a. Depending on the experimental setup, the imaged field of view will differ. In this example, source and the detector fibers are placed a certain distance ($\rho$) from each other and are coupled to the laser and detector. The imaged field-of-view (imaged over $i \times j$ pixels includes the fiber core which in later steps will be used to calculate $\kappa ^2$ over a specified region of interest ($w_z$). Sub-figure b illustrates Step 1 of the simulations. In this step, the rate at which the speckles decorrelate, $\tau _c$, is determined from the correlation diffusion equation (CDE). Using this value of $\tau _c$, consecutive frames of correlated speckles are simulated so that their electric-field autocorrelation decays with $\tau _c$. The intensity of these simulations are in arbitrary units, and independent of exposure time, $T$. Instead they represent speckles measured during a finite time-bin width, $t_{frame}$, on the $g_1$ curve. In order to simulate several values of $\rho$, the process illustrated in b can be repeated several times to simulate the $\rho$ dependent change in $\tau _c$. In Step 2 (sub-figure c), the arbitrary units of the simulated frames is scaled to represent realistic values of photon current rate, $\Phi$, in units of photons/second. In Step 3 (sub-figure d), an exposure time is introduced to the simulations by summing over frames. This process additionally converts the units of the simulations from photons/s to photons. Various values of $T$ can be simulated from the same set of simulated frames of Step 1. In this case, the simulation of two values of exposure time, $T_X$ and $T_Y$, is shown. Multiplying the summed frames in units of photons by the quantum efficiency (QE) of the camera converts the units of the simulations to electrons (e$^{-}$). In Step 4 (sub-figure e), the detector effects are simulated by altering the simulated intensity statistics according to the specifications of real detectors. In Step 5 (sub-figures f and g), $n$ speckles are sampled over an area, $w_z$ or over pixels of several repetitions of simulations to estimate a value of $\kappa ^2$. The yellow dots represent $\kappa ^2$ simulated for the $\tau _c$ and therefore $\rho$ simulated in Step 1. The two values of $T$ simulated in Step 3 are also shown. In the final step (Step 6, sub-figure h), the discrepancies in the exact form of the speckle autocorrelation decay between the solution for the CDE for a semi-infinite medium and the developed model is corrected for (NB: The scale of sub-figure h was modified from sub-figures f and g to emphasize the difference between the copula model and CDE.).

Download Full Size | PDF

2.1 Simulated experimental setup

Let us begin by detailing the canonical experimental setup that is being simulated. The exact details of the desired experimental setup to simulate may differ, however, the simulations are largely independent of these details. A visual representation of a possible setup is shown in Fig. 1(a). Here, the light is delivered through an optical fiber, and detected with a separate fiber coupled to a camera. The core of the fiber is imaged with appropriate optics and all the pixels within that region-of-interest (ROI) correspond to one value of source-detector separation, $\rho$. In a free-space system, the pixels in the imaged field of view could correspond to different values of $\rho$.

We assume that a coherent light source of wavelength $\lambda$ is utilized. The photons, once in the medium, undergo absorption and scattering events. The probability per unit length the photons are absorbed is estimated by the absorption coefficient ($\mu _a(\lambda )$). The reduced scattering coefficient ($\mu _s^{'}(\lambda )$) is used to estimate the total length which after a few scattering events leads to the randomization of the photon direction. In other words, after a photon traverses a distance few times the $1/\mu _s^{'}$, the light can be considered diffuse [42]. This diffuse light is measured at a distance $\rho$ away from the source. As a rule-of-thumb, $\rho$ is related to the mean probed depth by the measured light so that in order to measure deeper tissue, canonical experiments utilize longer $\rho$.

If the light source is of sufficiently narrow bandwidth (long coherence length) [43], then the so-called “diffuse laser speckles” and their statistical fluctuations can be observed. The electric-field ($g_1$) or the intensity ($g_2$) autocorrelation of the detected speckles are functions of parameters related to the experimental setup (e.g. $\rho$ and $\lambda$) and the properties of the measured medium including $\mu _a$, $\mu _s^{'}$, the ratio of the moving scatterers to the static ones ($\alpha$) and the mean-squared displacement of the scatters ($\Delta r^2$). For most experiments, the “effective” particle/scatterer diffusion coefficient weighted by $\alpha$ ($\alpha Db$) is measured as a “blood flow index” (BFI). For further details see Refs. [7,10,44]. The decorrelation time, $\tau _c$ (normally defined as the time $g_1$ decays to $1/e$ [20]) was defined for the purpose of these simulations as the time at which $g_1$ decayed to 0.5 and is also a function of these parameters.

2.2 Speckle statistics detected by a two dimensional detector array

We have simulated $\kappa ^2$ for tissue with specific optical properties and blood flow by simulating consecutive frames of correlated speckles which simulate their electric field autocorrelation with a decorrelation time, $\tau _c$, defined by the solution of the CDE for a semi-infinite medium [10]. The methodology presented is independent of this solution and other solutions (layered, heterogeneous, numerical) of the CDE could be utilized. For clarity, electric-field autocorrelation curves following the solution of the CDE will be referred to as $\hat {g_{1}}$, while the simulated electric-field autocorrelation curves are referred to as $\overline {g_{1}}$. While the two are similar, there are slight differences which are discussed below. Furthermore, the theoretical value of $\kappa ^2$ derived from the CDE will be referred to as $\hat {\kappa ^2}$ while the simulated values will be referred to as $\overline {\kappa ^2}$.

In the first step of the simulation pipeline (Figure 1(b)), $\tau _c$ is derived from $\hat {g_1}$. The derived value of $\tau _c$ was used to simulate frames of individual speckles by modifying the copula method developed in Ref. [45]. This method simulates consecutive two dimensional matrices of numbers that are correlated to each other by using a mathematical copula. Furthermore, the statistical profile of each matrix reflects the probability distribution of speckle intensity. Therefore, each individual matrix can be considered as a camera frame acquired in a speckle contrast experiment. These matrices are referred to as “frames” ($f$) simulating pixel coordinates $i,j$ while imaging speckles with diameter, Ø. Ø behaves as a scaling factor to put physical units for the pixel size since the speckle diameter is approximately equal to the wavelength of light being used. Therefore, choosing Ø to be equal to three pixels for a system modeling $\lambda = 785$ nm will scale the width of each pixel to be equal to approximately $262$ nm. The details of the simulation of Ø is described in the original copula method of Ref. [45].

The autocorrelation, $\overline {g_{1}}$, of the first frame, $f=f_1$ to the k^th frame, $f=f_k$ is given by

We note that the definition of $\tau _c = 0.5$ was an arbitrary decision. Depending on the definition of $\tau _c$ used, Eq. (2) can be updated to reflect the chosen definition. Each of the individual simulations of $\overline {g_{1}}$ consisting of $f=f_{N}$ frames of speckles patterns constitute an experiment, defined by $\epsilon$. $f_{N}$ frames are simulated to simulate the full $g_1$ curve to later simulate $\kappa ^2$ by integrating $\overline {g_{1}}$ from 0 to T [46]. This process together with notation is illustrated in Fig. 2. The basic method simulates $\beta$, an experimental parameter related to the coherence of the light source and the detection optics [47], equal to one. However $\beta$ can also be simulated for other values by following the method of Ref. [45].

 

Fig. 2. Illustration of how frames with a defined $\tau _c$ are simulated. First individual speckles are simulated on a grid of $i \times j$ pixels. These individual frames, $f$, are correlated to each other and their electric-field autocorrelation, $\overline {g_1}$, decay according to $\tau _c$ defined from semi-infinite theory (Figure 1). One full simulation of a theoretical $g_1$ curve ($\overline {g_1}$) consisting of $f_N$ frames corresponds to one experiment, $\epsilon$. This process is repeated several times resulting in several simulations of ${g_1}$.

Download Full Size | PDF

The simulations are simulated in arbitrary copula units. In addition, the frames are only dependent on $\rho$ and every simulated frame represents a point on the $\overline {g_{1}}$ curve with a finite time-bin width, $t_{frame}$. Since each frame has a defined $\rho$ and is simulated over an array $i\times {j}$, the complete notation is, ${}^{c}_{\sim }S(\rho )_{ijf}$. In this notation, the pre-superscript indicates the units of the simulated frame. In this case, $c$ refers to the arbitrary “copula” units. The pre-subscript, $\sim$, indicates that no effect of detector noise has been included in the simulated frame. The indices $i, j$ and $f$ refer to the pixel and frame.

2.3 Scaling detected photon intensity

The simulation of ${}^{c}_{\sim }S(\rho )_{ijf}$ is only a simulation of the decorrelation of speckles and is in arbitrary units. Therefore, in order to convert ${}^{c}_{\sim }S(\rho )_{ijf}$ to physical units and convey a realistic decay in detected photon intensity with $\rho$, the arbitrary copula units must be scaled to a realistic value (Fig. 1(c)). This is done by defining the spatial decay of light intensity theoretically or experimentally. According to the photon diffusion theory, in a semi-infinite geometry, the measured photon current rate, $\Phi (\rho )$, in units of photons/second, decreases with $\rho$ as:

Alternatively, experimental values of $\Phi (\rho )$ can be used to simulate the photon current rate at the detector. In this case, the average measured photons per second at specified values of $\rho$ (divided by the quantum efficieny of the specified detector) can be used to approximate the photon current rate.

Once $\Phi (\rho )$ has been established, whether theoretically or experimentally, the simulated frames are scaled using $\Phi (\rho )$ to convert them to a physically meaningful unit of photons/second, denoted as ${}^{ps}_{\sim }S(\rho )_{ijf}$. This is evaluated through the normalization of ${}^{c}_{\sim }S(\rho )_{ijf}$ with its mean over simulated frames, $\mu \left ( {}^{c}_{\sim }S(\rho )_{ij}\right )_f$:

2.4 Introducing exposure time to the simulated frames

The next step (Fig. 1(d)) requires converting the frames of equal frame widths, $t_{frame}$, to frames with an exposure time, $T_x$. These frames are denoted as ${}^{p}_{\sim }S(\rho, T)_{ijf}$ and are in units of photons. This is done by adding $N = T_x /t_{frame}$ consecutive frames:

Note that with the introduction of exposure time, the simulated frames drop their indexing of $f$. Furthermore, we also note that in order to properly introduce exposure time, it is important to sum simulated instantaneous frames with a sufficient time resolution. In other words, the integration of $g_1$ over $\tau$ (Ref. [46]), must be done with a small enough step-size. For example, for the case of $\beta = 1$, adding only the first frame with the last frame results in a simulated value of $\kappa ^2$ close to 0.5 independent of the value of $\tau _c$.

Finally, the simulated frames are converted from photons to electrons:

Table 1 summarizes the introduced notation to refer to the simulated frames.

Table 1. Table of definitions of the simulated speckle patterns including conversion of units from arbitrary simulation units with no $T$ dependency to electron units with $T$ dependency. In the notation for the simulated frames, the pre-superscript indicates the units of the simulated speckle intensities while the pre-subscript, $\sim$, indicates that no noise has been added

View Table | View all tables in this article

2.5 Detector noise

The final step before using the simulations to calculate $\overline {\kappa ^2}$ is to simulate the effects of the main types of detector noise on the simulated frames previously described, namely: photon shot noise, dark signal non-uniformity (DSNU), dark current shot noise, and read-out noise [48,49]. This step is illustrated in Fig. 1(e). To simulate detector noise, the distribution of each of the types of noise is considered, and random numbers are generated following the distribution. The notation used to describe the generation of random numbers and their distributions is shown in Eq. (7)

$I_{Z_{ij}}$ is the random number generated representing a certain intensity (in e$^{-}$) at pixel $i,j$. $I_{Z_{ij}}$ originates from a distribution, $p_Z$, with a mean value of intensity, $\mu (I)$, and variance, $\sigma ^2(I)$.

Photon shot noise is a Poisson distributed noise source [48,50]. Using the notation in Eq. (7), the contribution of photon shot noise at each pixel $i,j$ is described as:

DSNU and dark current noise along with read-out noise are not directly applied to ${}^{e}_{s}S(\rho, T_x)_{ij}$, instead independent dark frames are simulated and then added to ${}{_s^e}S(\rho, T_x)_{ij}$.

DSNU is simulated by simulating individual pixels of logistically distributed random numbers [50]:

The dark shot noise, similar to the photon shot noise (Eq. (8)) is simulated by applying Poisson distributed random numbers [48] to each pixel simulated in Eq. (9):

Finally, read out noise is simulated by assuming that it is a normally distributed noise source [51]. Read out noise in CMOS cameras is added at each pixel and is independent of the dark noise and the detected signal. Therefore, the contribution of the read out signal at each pixel, $I_{r_{ij}}$, is simulated:

The total dark frame, $df$, is then given by

Putting everything together, the frames with shot noise, DSNU, dark shot noise, and read-out noise, ${}^{\;\;e}_{sdr}S(\rho, T_x)_{ij}$, are given by:

To generalize the notation, the pre-subscript $N$ indicates a general noise source. In other words, ${}^{e}_{N}S(\rho, T_x)$ is shorthand for speckle intensity frames in units of electrons with unspecified noise, $N$, added. $N$ can take values:

The definitions and notation for simulating detector noise is summarized in Table 2:

Table 2. Table of definitions of the noise sources that are included in the simulations along with their corresponding distributions. The notation $p_Z(z; \mu, \sigma ^2)$ is used to define random numbers, $z$, originating from a distribution, $p_Z$, with a mean value of, $\mu$, and variance, $\sigma ^2$.

View Table | View all tables in this article

2.6 Speckle contrast

The final steps of the simulation pipeline require the calculation of $\overline {\kappa ^2}$ using the frames that have been simulated. In the first step, $\overline {\kappa ^2}$ is directly calculated using the simulated frames. The calculation of $\overline {\kappa ^2}$, as in a real experimental setting, can be done temporally or spatially depending on how speckles are sampled. Independent of the domain in which $\overline {\kappa ^2}$ is simulated, it should be noted that since the speckle decorrelation was modelled as a single exponential (Eq. (1)), the physically more realistic semi-infinite model of the speckle decorrelation follows a double exponential model [10]. A correction was applied in order to simulate a model corrected value of $\overline {\kappa ^2}$ denoted as $\overline {\kappa ^2}'$. Previous work in developing a successful DCS noise model also applied a single exponential model in order to model noise [28,52]. Therefore, while the value of $\overline {\kappa ^2}$ will be affected by the model used for $\overline {g_{1}}$, the noise is well described using the simplified single exponential model. The definitions and notation related to $\overline {\kappa ^2}$ are summarized in Table 3. The following sections will describe their calculations.

Table 3. Table of definitions for $\kappa ^2$. Three different variations of $\kappa ^2$ are calculated: first $\kappa ^2$ calculated directly from the integration of the double exponential $g_1$ from CDE. This is $\hat {\kappa ^2}$. Secondly, $\kappa ^2$ calculated directly from the simulated frames whose $g_1$ ($\overline {g_1}$) follows a single exponential form. This is $\overline {\kappa ^2}$ and outlined in Section 2.7. Thirdly, the model differences due to the differences in $g_1$ is corrected. This is $\overline {\kappa ^2}'$ and is outlined in Section 2.8. Moreover, $\overline {\kappa ^2}$ and $\overline {\kappa ^2}'$ can be calculated either spatially or temporally.

View Table | View all tables in this article

2.7 Model uncorrected speckle contrast

So far the process for simulating the detection of speckle statistics on a 2D detector array and the detector properties (Fig. 1(b) to (e)) has been described. These steps can be repeated in order to simulate several experiments ($\epsilon$, Fig. 2) for several different values of $\tau _c$ and therefore $\rho$, for calculating $\overline {\kappa ^2}$ in the temporal domain over $w_t$, or for determining $\sigma ({{\overline {\kappa ^2}}})$.

The next step in the pipeline is to use these frames to calculate values of $\overline {\kappa ^2}$ (Fig. 1(f) and (g)). As mentioned previously, $\overline {\kappa ^2}$ can be measured spatially or temporally i.e. speckle statistics can be determined spatially by using an area, $w_z$, of pixels or temporally over the pixels in a set of experiments, $w_t$.

Spatial $\overline {\kappa ^2}$ is given by:

Similarly, temporal ${\overline {\kappa ^2}}$ is given by:

With ${}_{N} {\overline {\kappa ^2}}$ simulated, noise correction must be applied. To do this, the noise correction method outlined in [2] was used. Here we outline the correction for spatial ${}_{N}{\overline {\kappa ^2}}$, but the same principles apply for temporal measurements.

Briefly, in order to correct for the dark and read signal offset in ${}_{N}{\overline {\kappa ^2}}$, a new dark frame, $df_{corr}$, is simulated using Eq. (12). The new dark and read signal offset corrected speckles frames is given by:

After the dark frame offset is corrected, the additional variance due to shot ($\sigma ^2_{shot}$) and the dark frame (dark and read out noise, $\sigma ^2_{df}$) is corrected by subtracting their respective variances from the signal variance, $\sigma ^2_{signal} = \sigma ^2({}^{\;\;\;\:e}_{sd'r'}S(\rho, T_x)_{w_z})_\epsilon$.

Putting everything together, the shot, dark, and read noise corrected value of ${\overline {\kappa ^2}}$, i.e. ${{}_{s'd'r'}{\overline {\kappa ^2}}_{{w_z}\epsilon }}$, is given by:

Variations in the noise correction can also be simulated. For example, the shot noise only added frames, ${}{_s^{}}{\overline {\kappa ^2}}$, can be corrected in the following way:

2.8 Model corrected speckle contrast

In these simulations, two forms of the electric field autocorrelation function have been introduced: $\hat {g_{1}}$ and $\overline {g_{1}}$, and crucially the decorrelation of the latter was modeled from the decorrelation time of the former. However, the two are described by two different exponential functions meaning that the values of $\kappa ^2$ derived from the two will differ. In particular, $\hat {g_{1}}$ describes a measurement in a semi-infinite medium and a multi-scattering (diffuse) regime. Since $\hat {g_{1}}$ is a more realistic solution to the CDE, rather than working with ${\overline {\kappa ^2}}$ derived from $\overline {g_{1}}$, we introduce another variable, ${\overline {\kappa ^2}}'$, which is the model-corrected value of $\overline {\kappa ^2}$.

${\overline {\kappa ^2}}'$ is derived from both ${\overline {\kappa ^2}}$ and $\hat {\kappa ^2}$. ${\overline {\kappa ^2}}$ values are used to simulate the offset or bias ($\gamma$) in $\kappa ^2$ due to noise, as well as to simulate the expected variance of $\kappa ^2$ over $\epsilon$. The CDE solution of $\hat {\kappa ^2}$ is then used to scale the value of ${\overline {\kappa ^2}}'$ to the expected value of speckle contrast when measuring in a semi-infinite geometry.

The bias term, $\gamma$ is defined as:

Finally ${\overline {\kappa ^2}}'$ values are generated by generating normally distributed random numbers, $k$, with mean equal to $\hat {\kappa ^2}+ \gamma$ and variance equal to $\sigma ^2({}{_{N}^{}}{\overline {\kappa ^2}})_{\epsilon }$:

No model correction was performed for the correction of the variance of ${\overline {\kappa ^2}}'$, following the same reasoning as was used in the well- validated noise models for DCS [28,52]. These previously developed DCS models have shown that the simulated variance in $g_1$, despite being simulated using a single exponential form of $g_1$, is equivalent for measurement geometries where $g_1$ takes a double exponential form [28,53–55].

2.9 Using the simulations to evaluate system performance

A primary motivation for developing a speckle contrast model is to evaluate the performance of such systems. Performance of simulated systems has been evaluated by its accuracy and precision. In this context, accuracy refers to the percent error of ${}_{N}{\overline {\kappa ^2}}'$ from its CDE solution, $\hat {\kappa ^2}$, and was defined as $100 \times \frac {{\overline {\kappa ^2}}' - \hat {\kappa ^2}}{\hat {\kappa ^2}}$. Precision is a measure of how variable a repeated measurement is and has been evaluated by its coefficient of variation (CV) as a percentage defined as the ratio of standard deviation of repeated experiments of ${}_{N}{\overline {\kappa ^2}}'$ to its mean: $100 \times \frac {\sigma ({}_{N}{\overline {\kappa ^2}}')_{\epsilon }}{\mu ({}_{N}{\overline {\kappa ^2}}')_{\epsilon }}$. Maximum accuracy and maximum precision correspond to the minimum values in these metrics.

2.10 Experimental setup (A) to validate simulations

The speckle contrast noise model was validated by comparing experimental results to the simulated noise for a range of exposure times. A multi-mode fiber delivered light (785nm, Crystalaser, Reno NV, USA), onto a a homogeneous, steady-state phantom with lipid droplets undergoing Brownian diffusion as in Ref. [29] was prepared. The resulting speckle pattern was imaged onto an sCMOS camera (Orca Fusion-C14440-20UP, Hamamatsu Photonics K.K., Hamamatsu, Japan) using a multi-mode fiber (910 $\mu$m core, 0.22 NA) and objective lens (f = 11 mm). The value of $\beta$ was measured to be approximately 0.2, and Ø was adjusted to be approximately 4 pixels.

$\tau _c$ of the system was obtained by simultaneous recording $g_2$ of the system using a single mode fiber coupled to a standard DCS device. The detector fibers of both the SCOS system as well as the DCS system were placed at a distance $\rho = 0.8$ cm from the source. The performance of the simulations was compared to the experimental results by evaluating the standard deviations of ${}_{sdr}\kappa ^2$ of both over 100 experiments. In addition, the expected signal-to-noise-ratio (SNR) was also evaluated considering $\mu ({\overline {\kappa ^2}}')$ to be equal to the average value of ${}_{sdr}\kappa ^2$ over 100 experiments (Eq. (20)). SNR is defined as the ratio of the average value of the signal over the noise. The experimental values of ${}_{sdr}\kappa ^2$ was calculated over a horizontal row of 1032 pixels. The simulated SNR was defined as the ratio of the standard deviation of the experimentally obtained values of ${}_{sdr}\kappa ^2$ to the average value of ${}_{sdr}{\overline {\kappa ^2}}'$ over 100 simulated experiments, $\epsilon$, calculated over 1032 simulated pixels.

2.11 Experimental setup (B) to optimize and design a speckle contrast system

The speckle contrast noise model was further used to design a speckle contrast system and define the required detected electron count rate (e$^{-}$/pixel/second) in order to accurately measure blood flow in the adult human brain. An sCMOS camera by Basler (daA1920-160um, Basler AG, Ahrensburg, Germany) was considered and simulated due its lightweight (15 g), compact size (19.9 mm x 29.3 mm x 29 mm) and cheap price (<300€). Measurements were chosen to be taken at $\rho$ of 2.5 cm and $T$ of 5 ms.

The required detected electron count rate to accurately measure $\kappa ^2$ was determined by attenuating a 785 nm laser (Crystalaser, Reno NV, USA) on a liquid using a fiber attenuator (OZ Optics, Ottawa Ontario, Canada). The diffuse light was imaged onto the camera using an 800 $\mu$m core multi-mode fiber (0.22 NA). The imaged speckles had a size of Ø = 5 pixels. The value of $\beta$ of the system was previously determined to be approximately 0.2. Speckle contrast data was acquired over 600 frames, and data was analyzed using an ROI of approximately 1100 pixels.

As in the setup (A) to validate the simulations, $\tau _c$ of the simulations was obtained from $g_2$ recorded using a standard DCS device. In order to approximate the required detected electron count-rate (e$^{-}$/pixel/second), a homogeneous (i.e., without layers) liquid phantom [29] was prepared to have optical properties of $\mu _a =$ 0.1 cm$^{-1}$ and $\mu _s^{'} =$ 10 cm$^{-1}$. The true value of $\kappa ^2$ was considered to be the value of $\kappa ^2$ measured with the highest detected intensity count rate, $I_{max}$. Percent error of $\kappa ^2$ as a function of the attenuated detected intensity count rates, $I_{att}$, was therefore calculated as: $100 \times \frac {\kappa ^2(I_{att}) - \kappa ^2(I_{max})}{\kappa ^2(I_{max})}$.

3. Results

3.1 Verification with experimental data

The results of the simulation model were compared to experimental data of an Orca Fusion camera using the experimental set-up in Section 2.10. Details of the camera parameters are summarized in Table 4. The simulations used $\tau _c$ obtained from the $g_1$ curve recorded using DCS (Figure 3(a)). $\beta$ was simulated to be 0.2 and Ø was set to 4 pixels to agree with the values of $\beta$ and Ø of the experimental data. Both experimental and simulation results were obtained for exposure times ranging between 0.1 ms and 5 ms in order to cover a range of detected electron intensities. It was ensured that the average value of the simulated detected electron intensity matched the experimental data (Figure 3(b)). The resulting experimental and simulated standard deviation of ${}^{e}_{sdr}\kappa ^2$ is shown in Fig. 3(c). The calculated signal to noise ratio of $\kappa ^2$ in Fig. 3(d), shows good agreement of the simulations with the experimental results.

 

Fig. 3. Comparison of the developed speckle contrast noise simulation model with experimental values. The number of experiments as well as the number of speckles used to obtain $\kappa ^2$ were the same for experiments and simulations. a) Experimental $g_1$ curves measured with a DCS system from which $\tau _c$ used in the simulations was determined (red). b) Average detected electrons over 1032 pixels and 100 experiments (black) and 100 simulations over 1000 pixels (grey). c) The standard deviation in ${}^{e}_{sdr}\kappa ^2$ calculated by simulation (grey) and the experimental results (black). d) SNR from experiment (black) and simulation (grey).

Download Full Size | PDF

Table 4. Simulation parameters used to verify simulations with experimental data acquired using an sCMOS camera (Orca Fusion-C14440-20UP, Hamamatsu Photonics K.K.)

View Table | View all tables in this article

3.2 Simulation study

Using the simulation pipeline described, we simulate speckle patterns with realistic detector noise. All simulations considered hardware consisting of a 785 nm unpolarized laser ($\beta = 0.5$) and a 100$\times$100 pixel array detector with noise properties derived from an Orca Flash4.0 v3 CMOS camera [56]. Since the variance of read-out noise is typically not defined in specification sheets, an online simulation tool was used to approximate the value of $\sigma ^2(I_r)$ [57]. Tissue with optical properties listed in Table 5 were simulated. These values were chosen as they are roughly the expected values when measuring in human tissue. $\overline {g_{1}}$ was simulated for $\rho$ ranging from 0.5 to 4.5 cm for $T_{max} = 5$ ms. Ø was chosen to equal three pixels in order to meet the requirements of the Nyquist criteria [38,58]. The details of the parameters used in the simulation are summarized in the table below:

Table 5. Parameters that were used to simulate synthetic speckles. Optical properties were chosen to mimic biological tissue, and detector parameters are based off of the properties of the Orca Flash4.0 v3 CMOS camera by Hamamatsu K.K.

View Table | View all tables in this article

3.3 Part I: simulating $\overline {\kappa ^2}$

The simulated values of the decorrelation time, $\tau _c$, as a function of source-detector separation, $\rho$, is shown in Fig. 4(a). As expected from theory, the speckle autocorrelation decays faster with increasing $\rho$ [10], confirming that the modified copula method for simulating decorrelating speckle intensity replicates the expected dynamics from theory. In Fig. 4(b), $\overline {\kappa ^2}(\rho )$ calculated by integrating the simulated speckle electric field decorrelation curves, $\overline {g_{1}}$ (Eq. (1)) for three different exposure times is shown. As expected from theory, $\overline {\kappa ^2}$ decreases with increasing $\rho$ and increasing $T$.

 

Fig. 4. a) simulated values of $\tau _c$ in ms. A clear decrease in $\tau _c$ with increasing $\rho$ is seen. b) ${\overline {\kappa ^2}}$ at three different exposure times calculated from integrating the autocorrelation, $\overline {g_{1}}$, of the simulated speckles.

Download Full Size | PDF

The simulated detected number of electrons (${}{_{\sim }^{e}S(\rho, T)_{ij}}$) for different $\rho$ at two different $T$ for all 100 simulated experiments are shown in Fig. 5(a) and (d). Including detector effects in the simulations results in deviations of the average value and variance from the ideal detected electron intensity value. This effect is $\rho$ and $T$ dependent. For all values of $\rho$ and $T$, the average value of the electron intensity does not deviate from the ideal case when only shot noise is simulated (N: $s$). However, in the regime of lower detected electron counts originating from speckle signal, i.e. at longer $\rho$ and shorter $T$, there is an increased variance in the shot noise included detected electron intensity. Furthermore, at short $T$, it is seen that the addition of a dark frame (N: $sdr$) visibly leads to a deviation in the average value of the detected electron intensity at $\rho = 2$ cm, while the same deviation for higher $T$ is not observed until approximately $\rho = 4$ cm. This is explained by the properties of the camera that were simulated. In this case, the dark current, a $T$ dependent signal, was significantly smaller than the read out signal, a $T$ independent signal, for the exposure times shown ($\mu (I_d) = 6\times 10^{-6} e^{-}$ and $\mu (I_d) = 3\times 10^{-4} e^{-}$ for $T = 0.1$ ms and $T = 5$ ms respectively, compared to $\mu (I_r) = 2.5 e^{-}$). Therefore, while dark noise is a $T$ dependent noise source, the effect of adding a dark frame appears more significant at shorter $T$ due to the high read-out signal relative to the speckle signal. Subtracting a dark frame (N: $sd'r'$) corrects this deviation. However a dark frame subtraction does not correct the increase in variance of the detected signal due to shot, dark, and read-out noise terms.

 

Fig. 5. Simulation of $\overline {\kappa ^2}$ from the frames of synthetic speckles. a, d) $\Phi (\rho )$ for two different exposure times ($T = 0.1$ ms and $T = 5.0$ ms on the top and bottom rows respectively) for when no noise source are added are shown as well as for when noise sources are added and when a dark frame is subtracted. b, e), the values of ${}{_N^{}}{\overline {\kappa ^2}}$ for all 100 simulated experiments. c, f) In order to correct for differences in theory of $g_1$ between the double exponential form of the semi-infinite model from CDE and the single exponential copula model, a bias term $\gamma$ is calculated (Eq. (19)). These are shown for different variations of added noise, N, at the two simulated exposure times.

Download Full Size | PDF

These observations are carried through to Fig. 5(b) and (e) where the values of $\overline {\kappa ^2}$ are plotted. At shorter $\rho$ and for both values of $T$, simulation of detector effects show very little deviation from the ideal, no detector noise added case. However, with increasing $\rho$, there is a noticeable deviation, as expected from experiments [2]. In the case of addition of shot, dark, and read-out noise (N: $sdr$), it is seen that for $T = 0.1$ ms (Figure 5(b)), ${}_{sdr}{\overline {\kappa ^2}}$ begins to deviate from the ideal case, at approximately $\rho$=2.0 cm. At $T = 5.0$ ms (Figure 5(e)), ${}_{sdr}{\overline {\kappa ^2}}$ begins to deviate from the ideal case from approximately $\rho$=1.5 cm. Correcting for detector effects by applying a dark frame subtraction and correcting for shot, dark, and read-out noises (N: $s'd'r'$) results in a larger range of $\rho$ for which $\overline {\kappa ^2}$ agrees with the ideal case for $T$=5.0 ms, to about $\rho$=3 cm. However, the same correction does not obviously perform as well for $T$=0.1 ms (Figure 5(b)), with detector effects correction (N: $s'd'r'$) apparently performing worse than the uncorrected case (N: $sdr$). This last observation should not be interpreted as a failure in the correction of noise, rather it is a reflection of the origin of the electron signal in this regime. Referring back to the plot of the detected intensity (Figure 5(a)), at $T$=0.1 ms, the majority of the detected electron signal after $\rho$=2 cm originate from the detector rather than from speckles. Therefore, without applying corrections, any value of $\kappa ^2$ in this regime is not a reflection of speckle contrast, rather reflects a “detector signal" contrast.

The bias term, $\gamma$ (Eq. (19)), is shown in Fig. 5(c) and (f) and reflects the offset of ${}_{N}{\overline {\kappa ^2}}$ from the no noise added case, ${}_{\sim }{\overline {\kappa ^2}}$. These were used to calculate the average theory corrected value of $\kappa ^2$ with simulated detector effects (${}_{N}{\overline {\kappa ^2}}'$). For the remaining results, only the case of N = $s'd'r'$ will be considered as this is the case of most interest in any experiment. The theory corrected values of $\kappa ^2$ are shown in Fig. 6(a) and (d).

 

Fig. 6. a, d) Simulation of theory corrected values of speckle contrast, ${}_{\sim }{\overline {\kappa ^2}}'$. b, e) Accuracy (percent error) of ${}_{\sim }{\overline {\kappa ^2}}'$. c, f) Precision (coefficient of variation) of ${}_{\sim }{\overline {\kappa ^2}}'$.

Download Full Size | PDF

Theory corrected values of speckle contrast, ${}{_N^{}}{\overline {\kappa ^2}}'$, were calculated from Eq. (20). The final averaged value of the simulated 500 normally distributed random values of ${}{_N^{}}{\overline {\kappa ^2}}'$ for $T=0.1$ ms and $T=5$ ms are plotted in Fig. 6(a) and (d). Error bars reflect the standard deviation. The accuracy of ${}{_N^{}}{\overline {\kappa ^2}}'$ is shown in Fig. 6(b) and (e), reflected as the percent error. The percent error increases (accuracy decreases) with increasing $\rho$ reaching $5{\% }$ at approximately 1.8 cm for short $T$ (Fig. 6(b)) and $2.5$ cm for long $T$ ((Fig. 6(e)). Similarly, the precision of ${}{_N^{}}{\overline {\kappa ^2}}'$, represented as the coefficient of variation (CV) also decreases (CV increases) with increasing $\rho$ (Fig. 6(c) and (d) for $T=0.1$ and $T=5.0$ ms respectively).

3.4 Part II: Using the simulations to study precision and accuracy

As seen in the previous section, effects of detector noise lead to decreases in accuracy of $\overline {\kappa ^2}'$ particularly in the regimes of long $\rho$ and short $T$. In the next part of this analysis, the simulations are used to understand how various parameters can be changed in order to increase the usable range of $\rho$ and $T$ considering both precision and accuracy. In order to quantify the requirements of a SCOS or SCOT system, it is assumed that the required accuracy is within a 5% error and precision within a 10% coefficient of variation (CV) at $\rho$=4 cm and $T$=5 ms. These values were chosen for deep tissue measurement: $\rho$=4 cm corresponds to an approximate measurement depth of $2$ cm. Although $\rho$=2.5 cm is considered sufficient for measuring the cortical surface going to further distances offers greath depth sensitivity and distances of between 3.0 - 4.0 cm have been used for tomographic reconstruction of human functional activation [59,60]. $T$=5 ms was chosen in order to be able to sample at fast enough acquisition rates while also maximizing the number of detected photons (Figure 5(d)).

In speckle contrast optical tomography (SCOT) or speckle contrast diffuse correlation tomography (scDCT) [16,17], several source and detector positions are used in order to reconstruct a three dimensional image of blood flow. In a system incorporating nine source positions as in [61], using $T$=5 ms, this will correspond to a full acquisition rate of 22.2 Hz for $\kappa ^2$ measured at each source position. Furthermore, 5% accuracy and 10% precision have been chosen as our targets since a 10% blood flow change corresponds to approximately 10% change in $\kappa ^2$. A 10% change in flow is similar to what is measured in functional studies [23].

It is known that a contributing factor to the precision of $\kappa ^2$ is the number of speckles used to determine $\mu$ and $\sigma ^2$ [34,38]. In the previous simulations of $\overline {\kappa ^2}$, $w_z = 100 \times 100$ pixels corresponding to the sampling of 1100 independent speckles. In Fig. 7, $w_z$ was changed to simulate the effects of the number of independently sampled speckles on the CV of $\overline {\kappa ^2}'$.

 

Fig. 7. The effect of the number sampled speckles on the measured precision of ${}{_{s'd'r'}^{}}\overline {\kappa ^2}'$ at three values of $\rho$, and $T$ = 5 ms. Increasing the number of sampled speckles results in a decrease in the CV of ${}{_{s'd'r'}^{}}\overline {\kappa ^2}'$.

Download Full Size | PDF

As expected in Fig. 7, increasing the number of speckles used to calculate $\kappa ^2$ results in an increase in the precision of $\kappa ^2$. The decay in CV with increasing speckle number follows a square root dependency, in accordance to the theory [34]. Therefore, if the objective is to measure $\kappa ^2$ with 10% precision at $\rho$=4 cm and $T$=5 ms, $w_z$ must be increased from 100 x 100 to approximately 170 x 170 pixels corresponding to approximately 3000 speckles (since Ø=3 pixels). Sampling more speckles can easily be implemented in a typical sCMOS camera with 2048$\times$2048 pixels by choosing a larger region of pixels.

As observed in Fig. 6(b) and (e), accuracy was seen to be higher at shorter $\rho$ and longer $T$, i.e. in the regime of high $\Phi$. Strategies for increasing the amount of detected light to achieve good accuracy while remaining within safety limits may include employing dual sources located equi-distance apart from the detected area of interest.

In addition to $\Phi (\rho )$, $\tau _c$, may also affect accuracy of $\kappa ^2$. In order to study the effect of $\tau _c$ on accuracy in $\kappa ^2$, the simulations were repeated fixing $\Phi (\rho )$ to be constant over all values of simulated $\rho$.

In Fig. 8, the percent error in ${}{_{s'd'r'}^{}}\overline {\kappa ^2}'$ as a function of the number of detected electrons shows that measurement accuracy is dependent on $\rho$, and by extension, $\tau _c$. For the simulated camera, measurements with longer $\rho$ (shorter $\tau _c$) require less detected electrons to achieve the same accuracy in $\kappa ^2$.

 

Fig. 8. Accuracy of ${}{_{s'd'r'}^{}}\overline {\kappa ^2}'$ for two different values of $\rho$ with identical values of $\Phi$ ($T = 1 ms$). Higher accuracy was found for greater $\rho$.

Download Full Size | PDF

The results of Fig. 8 indicate that besides camera properties and detected light intensity, variables affecting $\tau _c$ such as $\rho$ or the optical properties of the tissue can also affect the accuracy of $\kappa ^2$ measurements. We hypothesize that this result may be explained by the fact that the distribution of speckle intensity changes with $\tau _c$ (in particular the variance), and as a result there is a change in the sampling requirements for properly sampling the distribution. However, further studies should be performed to verify this hypothesis.

3.5 Using the simulations to design and optimize a system

In the previous sections we have verified the simulation pipeline by comparing the SNR measured experimentally with an Orca Fusion-C14440-20UP camera to the expectations from simulation. We have further demonstrated in detail (without experimental comparison) the entire simulation pipeline. Finally, in the following section we will demonstrate how these simulations can be used to design and optimize a speckle contrast system.

Speckles were simulated using the parameters specified in Table 6. These parameters were derived from the experimental results ($\tau _c$ and Ø), properties of the camera defined by the manufacturer, as well as data analysis ($w_z$). The resulting experimental and simulated percent error in $\kappa ^2$ for varying detected electron count rates is shown in Fig. 9.

 

Fig. 9. The effect of changing values of detected electron count rate on both the experimental and simulated values of percent error of $\kappa ^2$. The grey horizontal line marks 5% error.

Download Full Size | PDF

Table 6. Parameters that were used to simulate synthetic speckles based on experimental data taken using a Basler (daA1920-160um) CMOS camera on a liquid phantom.

View Table | View all tables in this article

The experimental and simulated results are in good agreement with each other and suggest that for the chosen detector, a minimum detected count rate on the order between 4 to 5$\times$10$^4$ e$^-$/pixel/second allows us to calculate $\kappa ^2$ with approximately 5${\% }$ error.

Using the derived acceptable minimum detected count rate as a guide in determining the accuracy of raw data signal, the same device was placed on a human subject’s forehead using a $\rho$ of 2.53 cm and $T$ of 5 ms. Data was acquired at a frame rate of 100 fps. A summary of the measurements is show in Fig. 10. The desired electron count rate was reached (around 4.3$\times$10$^4$ e$^-$/pixel/second, Fig. 10), and the resulting $1/\kappa ^2$ shows the expected pulsatile behavior for a measurement acquired at this frame rate (Fig. 10(a)). In order to confirm that the pulsatile behavior has physiological meaning, the fast Fourier transform (FFT) of the data has also been plotted (Fig. 10(c)). A distinct peak at 1.4 Hz is seen in the FFT corresponding to a heart rate of 84 bpm. This value matches the resting heart rate measured in this subject using a standard pulse oxymeter. The harmonics of the heart rate can also be seen.

 

Fig. 10. Summary of results from a SCOS measurement on an adult human forehead. a) $1/\kappa ^2$, a surrogate measure of blood flow, shows clear pulsatile signals. b) Average detected electron count rate lies in a range which allows us to accurately measure $\kappa ^2$. c) Fourier transform of the $\kappa ^2$ signal. A clear peak is found at 1.4 Hz corresponding to the heart rate of the subject (84 bpm).

Download Full Size | PDF

4. Discussion

A comprehensive model of speckle contrast signal for measurement of flow requires three main components: the simulation of speckles, their dynamics, and the detector effects on the measured signal. Individual 2D frames of speckles with the correct intensity distribution in these simulations were simulated following the method of Duncan et.al. [62]. The dynamics of the speckle intensity were simulated modifying the method of Ref. [45], where crucially the modification allowed for the characterization of $\tau _c$ to be specified according to speckle intensity decorrelation defined by the correlation diffusion equation [10]. While the exact form of the speckle decorrelation, $g_1$, differs in the simulations, general properties of the dynamics and their dependency on parameters such as $\rho$ and $\alpha Db$ could be simulated. The simplification of $g_1$ of a semi-infinite medium as a single exponential function has been seen to be accurate in noise models for DCS [28]. We note that this work-flow is applicable to beyond the photon diffusion regime and in the presence of static scatterers through the use of the appropriate model. Detector effects were simulated taking into account photon shot noise, dark current signal and noise, and read-out signal and noise. Our method for modeling speckle contrast can account for parameters such as the speckle to pixel size and $\beta$.

The developed model was validated experimentally, and was shown to accurately predict the SNR of $\kappa ^2$ measurements for a wide range of exposure times and detected intensities (Fig. 3). We noticed a greater discrepancy between simulations and experiments for smaller exposure times where the detected intensity is very small. We hypothesize that in this regime of low detected intensity, this discrepancy can be a result of having neglected some sources of detector noise such as quantization error. In addition, for regimes with noisier $\kappa ^2$, larger sample sizes (i.e. greater $\epsilon$) is required for accurately calculating $\sigma (\kappa ^2)$.

In the simulation study to evaluate system performance, we have shown that the simulations accurately represent experimentally observed behavior of $\kappa ^2$ in the regime of long $\rho$ and/or short $T$ where the speckle contrast signal increases above the theoretically expected values. Simulation of the noise correction method of Ref. [2] extends the region of $\rho$ and $T$ where the speckle contrast signal matches its theoretical value. However, depending on the amount of the contribution of the detector effects, the correction cannot account for all of the increased variance from these effects. Therefore, it is important when designing a speckle contrast system to consider the range of $\rho$ and $T$ where $\kappa ^2$ can be corrected. We have also shown the dependency of accuracy in speckle contrast signal on parameters including the number of detected photons, $\rho$, and $\tau _c$.

The accuracy and precision of $\kappa ^2$ developed in the simulation model not only reflects observed experimental behavior, but is also comparable to what has been described in the noise models of related techniques. In DCS, similar to what we have seen in speckle contrast, the SNR of the raw $g_1$ signal is dependent on the detected photon intensity and $\tau _c$. Since DCS uses correlators to measure $g_1$, the noise model for DCS also depends on the architecture of the correlator [28,63]. An emerging variation of DCS known as interferometric DCS, or iDCS, utilizes a heterodyne detection technique mixing the traditional DCS signal with a reference arm (i.e. the coherent source). This detection scheme results in greater values of $\tau _c$ compared to traditional DCS resulting in an increase in the SNR of the raw $g_1$ data as well as a decrease in the coefficient of variation of the retrieved blood flow values [15].

While in this analysis we have concentrated on the effects of detector noise in the regime of low detected photon counts corresponding to the typical observations in experiments, it is worth noting that high photon count rates that saturate the detector can also lead to decreases in accuracy as well as precision of the raw signal and in the derived blood flow values. In DCS, saturated detection leads to decreases in the experimentally measured $\beta$ resulting in inaccuracy of the retrieved blood flow [32]. Although not shown here, the same applies in measurements of speckle contrast as detector saturation will lead to inaccurate measurements of $\sigma ^2(I)$ and/or $\mu (I)$ and consequently $\kappa ^2$.

The copula method [62] has previously been used by Qiu et.al. [35] to study the effects of pixel sampling (sampling of $w_z$ and $w_t$) on $\kappa ^2$. In this work, a pseudo exposure time was considered. However since the decorrelation of the speckles were not reassigned in units of time, the simulations were not related to proper physiological or system properties. Thompson et.al [37] combined the method of simulating a single frame of speckles of Ref. [62] with small random phase changes for each consecutively simulated frame, making it very similar to the copula method of Ref. [62]. These simulations were used to study the effect of speckle to pixel size ratio in the measurement of $\kappa ^2$. However, like in Ref. [35], the simulations were not scaled to represent physiological properties and did not include any effects of detector noise.

The present study is complementary to the recent publication by Zilpelwar et.al. [40], with several notable differences. The model developed by Zilpelwar et.al. is based on a Monte-Carlo method simulating random particle (scatterer) motion. Their approach considers a single scattering regime, and is therefore strictly speaking is not applicable for SCOS which is a diffuse optical method considering a multi-scattering regime. Our approach does not simulate particle motion, rather we directly simulate the statistical properties of decorrelating speckle by generating correlated random numbers using the method of Duncan et.al. [45]. Both simulations are based on a single-exponential form of $g_1$. In the present work, we argue that while the exact value of $\kappa ^2$ is dependent on the approximations used to define $g_1$, the noise in $\kappa ^2$ is likely not affected due to previous observations in the development of a noise model for DCS [28]. In order to account for the difference in $\kappa ^2$ stemming from discrepancies in the approximation of $g_1$, in our simulations, we have included a method to correct for this difference. Furthermore, in the present work we were interested in deriving limits of accuracy and precision for an experimental scenario and therefore included a full noise corrected simulation of ${}_{\sim }{\overline {\kappa ^2}}'$ by simulating the expected dark frames of the individual specifications of each simulated camera. These details, multi-scattering regime in a semi-infinite medium, was not included in the model of Ref. [40].

Another recent publication by Murali et.al. also proposes a method for simulating decorrelating speckle statistics including the effect of detector noise [64]. In this work, unlike our proposed method or the method of Zilpelwar et.al. [40], the decorrelating speckle statistics are simulated for any form of $g_1$. This method can potentially offer better results for more complicated geometries or in experimental situations with significant contributions from static scatterers. However, for most typical SCOS or SCOT experimental scenarios, we expect that the proposed simple copula-based method is sufficient for modeling speckle contrast data. In the present work a detailed study comparing existing simulation methods such as those of Refs. [40,64] was not performed, however a direct comparison is warranted in the future.

We are not the first to adapt the work of Duncan et.al. [45,62] to study the behavior of $\kappa$. We note that this method is not only method in the literature for simulating decorrelating speckle patterns [65–68]. In the copula method of [45], spatial correlation is not preserved between frames. Song et.al propose another method for simulating frames correlated in the spatio-temporal domain [65]. The authors successfully simulated real speckle contrast data by creating correlation maps of data from a rat ear, however the authors note that the accuracy of replicating an image taken from real data depends greatly on the quality of the camera used to acquire the image. Sang et.al. utilized the method of Song et.al. [65] to further expand the method to include time integration effects of exposure time [69], however only one exposure time was simulated. Another method for modelling speckles is to model the summation of random phasors [66]. Postnov et.al. modified this technique in order to simulate the effects of the laser linewidth and camera noise on $\kappa ^2$ [67]. Finally, we also note that a separate method for simulating static speckle using a Monte-Carlo method has also been developed by Bar et.al. [70].

An interesting work by Song et.al. [71] derives the effect of camera quantization of intensity on speckle contrast from the probability density function of speckle intensity. Quantization of the speckle signal is something that was not considered in the current study and should be considered in future work.

5. Conclusion

In the present work we have introduced a method for simulating the formation and detection of dynamic speckle patterns. The main application that we have focused on was the design and characterization of a speckle a contrast system capable of measuring human adult cerebral blood flow non-invasively. To this end, the simulation method was validated on a dynamic liquid phantom, the details of speckle contrast signal as a function of $\rho$ and $T$ were studied, and finally a system designed for human cerebral blood flow was characterized and validated on an adult human subject.

Similar recent publications in the field highlight the need for methods for simulating speckle contrast signal and noise [40,64]. In contrast to these publications, the main contribution of the present work is the presentation of a full pipeline for the design and characterization of a SCOS system with clearly defined experimental requirements including the detection accuracy and precision of $\kappa ^2$. The simulation method has been shown to be useful when identifying the lower bounds of detected electron count-rate to achieve the desired accuracy and precision of speckle contrast signal. As speckle contrast signal is sensitive to detector noise effects at low detected electron count-rates, characterizing these limits is advisable when developing any speckle contrast system.