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

Laser speckle simulation tool based on stochastic differential equations for bio imaging applications

Open Access Open Access

Abstract

Laser speckle-based blood flow imaging is a well-accepted and widely used method for pre-clinical and clinical applications. Although it was introduced as a method to measure only superficial blood flow (< 1mm depth), several recently introduced variants resulted in measuring deep tissue blood flow (a few cm) as well. A means of simulating laser speckles is often necessary for the analysis and development of these imaging modalities, as evident from many such attempts towards developing simulation tools in the past. Such methods often employ Fourier transforms or statistical tools to simulate speckles with desired statistical properties. We present the first method to use a stochastic differential equation to generate laser speckles with a pre-determined probability density function and a temporal auto-correlation. The method allows the choice of apriori gamma distribution along with simple exponential or more complex temporal auto-correlation statistics for simulated speckles, making it suitable for different blood flow profiles. In contrast to the existing methods that often generate speckles associated with superficial flow, we simulate both superficial and diffuse speckles leading to applications in deep tissue blood flow imaging. In addition, we have also incorporated appropriate models for noise associated with the detectors to simulate realistic speckles. We have validated our model by comparing the simulated speckles with those obtained from in-vivo studies in mice and healthy human subject.

© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

Non-invasive imaging of microvascular blood flow using coherent laser source illumination and associated scattered speckle intensity [1,2] measurements have been studied over several decades with promising results. Conventionally, as in laser speckle contrast imaging (LSCI) [3] systems, the sample is illuminated uniformly and the speckle contrast of the measured speckles is used to quantify the surface blood flow up to a depth of 1mm. A spatio-temporal speckle contrast measured either at single exposure time [3] or multiple exposure time [4] of the camera are often employed for this purpose. With the advent of high-speed cameras and single-photon avalanche diode (SPAD) arrays, utilizing intensity autocorrelation instead of the time-averaged speckle contrast to quantify blood flow is gaining much attention [58].

On the other hand, diffuse correlation spectroscopy (DCS) uses a focused laser illumination to measure the deep tissue blood flow (order of a few cm) using spatially resolved diffuse speckle measurements. Speckle intensity measured at multiple source-detector (SD) separations are used to quantify blood flow, with the propagation model being the Correlation diffusion equation (CDE) [2,9]. Most often, a high-density measurement of speckles at several SD pairs is required either for better signal-to-noise ratio (SNR) [9] or to implement so-called called diffuse correlation tomography (DCT) [10]. However, in diffuse speckle contrast analysis (DSCA) [11], an inexpensive alternative to high-density autocorrelation measurements, speckle contrast from multi-speckle measurements from low frame rate cameras was used to measure deep tissue blood flow. Concurrently, speckle contrast optical spectroscopy (SCOS) [12] and its tomographic counterpart called speckle contrast optical tomography (SCOT) [13] were introduced with a better quantification of flow using spatially resolved speckle contrast measurements.

There have been other variants of the speckle based blood flow imaging such as interferometric DWS (i-DWS) [14], speckle contrast diffuse correlation tomography (sc-DCT) [15], multi-speckle DCS (M-DCS) [16,17], multi-exposure interferometric diffusing wave spectroscopy (mi-DWS) [18], Interferometric speckle visibility spectroscopy (ISVS) [19], Integrated Diffuse Speckle Contrast Spectroscopy (iDSCS) [20], holographic Fourier domain DCS [21], and inverse of the numerical integration of squared g1 (INISg1) [22]. In addition to the blood flow measurements, laser speckles have been employed to measure viscoelastic properties as in laser speckle rheology [23] and ultrasound-assisted DWS [24,25].

Although many such modalities use laser speckles to measure blood flow, a rigorous methodology for the simulation of laser speckles is lacking. It is not only essential for improving the quantification of blood flow by optimizing the source and detector characteristics but also for having a quantitative comparison among these modalities. Some of the existing methods used to generate laser speckles in the past are reported as follows. One common approach is using a fast Fourier transform of the phase matrix, which generates independent speckle patterns, but without any pre-determined correlation statistics [26]. In Ref [27], Duncan et al., proposed a method based on a copula algorithm for generating speckle sequences for a given auto-correlation model that utilizes the concepts of the quantile function and Fourier transform. However, this method was not utilized for generating spatially varying diffuse speckle patterns. The concept of coherent imaging principles was used to generate spatially varying speckle images using pre-defined correlation coefficient matrices in Ref [28]. A comprehensive review of the above-mentioned methods used for generating speckles can be found in chapter 3 of Ref [29]. Very recently, in Ref [30], James et al., have used a method based on the Karhunen–Loève expansion of the electric field to simulate time-integrated speckle patterns. All the above-mentioned work stress the need for laser speckle simulation tools for blood flow imaging applications.

A new approach for the simulation of speckle intensity based on stochastic differential equation (SDE) is introduced in this paper. A SDE, in simpler terms, is an ordinary differential equation with a stochastic source term, whose solution is a stochastic process. SDE has several applications in various fields that include, but not limited to, chemistry, finance, epidemiology, mechanics and microelectronics [31,32]. Here, for the first time, we apply the SDE to simulate laser speckles for blood flow imaging. The laser speckles originating from a turbid medium are characterized by probability density function (PDF) which is either negative exponential or gamma and with a simple exponential or complex temporal autocorrelation [33]. The coefficients of the SDE are determined in such a way that the above-mentioned statistics of the speckles can be incorporated into the simulation of speckles. In addition, spatially varying diffuse speckles are also simulated by appropriately modifying the coefficients of SDE using the solution to photon diffusion equation with a potential for application in deep-tissue blood flow imaging. In order to make this into a comprehensive laser speckle simulation tool, we have also applied realistic noise models to the simulated speckle pattern based on noise associated with the common speckle measuring detectors like array detectors (cameras) and avalanche photodiodes (APD).

We have compared the simulated speckles using SDE with the $in-vivo$ experimental speckles generated from the cerebral blood flow in mice and healthy adult humans. The presented speckle simulation tool incorporating both conventional and diffuse speckles with the flexibility to incorporate pre-determined PDF and complex temporal autocorrelations could guide researchers in choosing appropriate models and techniques to achieve better blood flow quantification.

2. Simulation model for speckle intensity using SDE

2.1 Theory

Consider a SDE which consists of an ordinary differential equation (ODE) along with a stochastic term [3436], given as

$$\begin{array}{r} dx(t)=a(x(t),t)dt + b(x(t),t)dW(t); \\ x(0)=x_0. \end{array}$$

Here a and b are called the drift and diffusion terms, respectively, and W is the Wiener process (or Brownian motion) such that $dW(t) \approx \sqrt {\delta t} N(0,1)$, where N denotes a Gaussian random variable with zero mean and unit variance. It is well known that functions $a$ and $b$ satisfy the Fokker-Planck equation [31], which is given as

$$\frac{\partial (ap)}{\partial x} = \frac{1}{2}\frac{\partial^2 (b^2 p)}{\partial x^2}.$$

Here $p$, PDF of random variable $x$, is assumed to be independent of time, as it is a stationary random process. Note that $a$ and $b$ implicitly depend on time. For $x(t)$ to represent speckle intensity, it should have a negative exponential PDF, i.e., $p(x)=\frac {1}{\mu } exp(-x/\mu )$, where the mean and the standard deviation are equal to $\mu$ [33]. Additionally, the auto-correlation of the speckle intensity obeys exponential decay [1,37]. The usual procedure to enforce this in the solution of Eq. (1) is to choose the parameters $a$ and $b$ using Eq. (2) appropriately. It can be seen that the terms $a(x(t)) = -\alpha (x-\mu )$ will ensure a correlation decay of the form $e^{-\alpha \tau }$ , which results in the equation $dx(t) = -\alpha (x - \mu )dt + b(x(t))dW(t)$ [34,35]. This predetermined $a$ is now substituted into Eq. (2) along with the exponential PDF to deduce $b$ as follows:

$$\begin{array}{r} b^2(x(t))= \frac{2}{p(x)}\int_0^x a(x(s)) p(s) ds \\ = \frac{2}{p(x)}\int_0^x -\alpha(x(s)-\mu) p(s) ds \\ = 2\alpha\mu x. \end{array}$$
The resulting SDE (with a change of notation $I \equiv x$), is given by
$$dI(t)={-}\alpha (I(t)-\mu) dt + \sqrt{2\alpha\mu I(t)} dW(t),$$
whose solution gives the speckle intensity with a prior PDF and auto-correlation. Note that the above SDE as in Eq. (3) is known as Cox-Ingersoll-Ross (CIR) model ($dX(t)=-\alpha (X-\mu )dt + \sigma \sqrt {X(t)} dW$) in mathematical finance, where $X(t)$ represent the interest rate. Clearly for our case, $I(t)$ has to be positive, which is ensured by Feller’s condition, given as $\frac {2\alpha \mu }{\sigma ^2} \geq 1$. In our case, $\frac {2\alpha \mu }{\sigma ^2}=\frac {2\alpha \mu }{2\alpha \mu }=1$ and hence $I(t) > 0$ for all $t$. The resulting SDE in intensity is solved to obtain the speckle intensity pattern.

2.2 Computational implementation

Although there are several numerical methods to solve Eq. (3), we chose to use the well known Milstein scheme [32,34,35,38,39]. The Eq. (3) solved by Milstien scheme is given by

$$I_{n+1} = I_n + a(I_n) \Delta t + b(I_n)\Delta W_n + 0.5 b(I_n)b^{'}(I_n)[(\Delta W_n)^2 - \Delta t],$$
with $I_0=\mu$ and $n$ is the index of time step.

The SDE algorithm used to generate speckles, $I(t)$, is shown as a pseudocode below in Algorithm 1.

Tables Icon

Algorithm 1. SDE - Milstein scheme

$*-$Input parameters; $\mu$ is mean intensity based on photon diffusion equation and $\tilde {\alpha }$ is from the auto-correlation function, based on the flow.

$\#-$randn - MATLAB command - generate random numbers

$*\#-$gives speckles with mean $\mu$ and correlation decay $1/\tilde {\alpha }$.

The autocorrelation $g_2$ was calculated using the standard MATLAB function $xcorr$ (with an unbiased estimate of the cross-correlation), wherein $g_2$ is defined as $g_2(\tau ) =\frac {(<I(t)I(t+\tau )>)}{((N-\mid m \mid )<I(t) >^2)}$, $N$ is the number of samples, and $m$ is the index of the lag $\tau$. This unbiased estimate of $xcorr$ is essential for obtaining a proper estimate of intensity autocorrelation, while biased and other forms of autocorrelation results in distorted or undesired autocorrelation for longer time integration. The effect of using biased/ unscaled autocorrelation is shown in Fig S3 of the supplemental document, wherein the biased auto-correlation does not follow the desired auto-correlation at a longer delay time, while unbiased auto-correlation follows the desired auto-correlation. Here $g_2^N=\frac {(<I(t)I(t+\tau )>)}{max(<I(t)I(t+\tau )>)}$ and $g_2^N (Biased/unscaled) =\frac {(g_2-1)}{\beta }$. This is due to the insufficient number of samples leading to an inherent bias at larger values of $\tau$, resulting in distortion of autocorrelation.

2.3 Simulation parameters

Depending on the application, we can take auto-correlation to be independent of spatial coordinates $r$ (in a sample with uniform illumination as in LSCI/ MESI) or dependent on $r$ as in diffuse speckles (DCS, DSCA, SCOS, sc-DCT, M-DCS and i-DWS). We have used $t_{min}=10^{-6} s,t_{max}=1s$ and $N=10^7$ samples to generate the speckle intensity.

2.3.1 Surface blood flow

For surface blood flow, we have normalized intensity auto-correlation $g_2(\tau )=1+\beta |g_1^2|$, where $g_1=e^{-t/\tau _c}$ and $\tau _c$ is the characteristic decay time (based on Siegert’s relation [2,40], $\beta$ is a constant depending on the collection optics of the experiment). In order to demonstrate the capability of the proposed method in simulating complicated autocorrelation models (henceforth referred to as complex autocorrelations) other than the simple exponential based one, we consider the autocorrelation model proposed in Ref [6] given by ,

$$\begin{array}{r} g_2(\tau)=1+\beta \rho^2 (d\vert g_1^{(n=X)}(\tau)\vert+(1-d) \vert g_1^{(n=1)}(\tau)\vert)^2 + \sqrt{\beta} (1-\rho)\rho \\ (d \vert g_1^{(n=X)}(\tau)\vert +(1-d)\vert g_1^{(n=1)} (\tau)\vert )+C \end{array}$$

Here $\beta$ is the instrumentation constant, $X$ and $d$ dictate the dynamic scattering regime, namely, multiple unordered (for $X=0.5$), single unordered or multiple ordered regime (for $X=1$) and single ordered regime (for $X=2$) and $\rho$ is the fraction of dynamic scatterers. Three representative tissue types were used for the simulations: large vessels, small vessels and parenchymas. The values of the above parameters for modelling were chosen based on Ref [6]. For the large vessels, $\tau _c=50 \mu s$ was used along with $\rho =0.98, d=0$ and $X=1$. For small vessels, we use $\tau _c=250 \mu s$, $\rho =0.98, d=1$ and $X=2$ and for the parenchyma $\tau _c=250 \mu s, \rho =0.8$, $d=1$ and $X=0.5$ was used.

2.3.2 Simulation of complex auto-correlation models

Any complex exponential function (such as $e^{k\sqrt {\tau }}$ or $e^{k \tau ^2}$ or its combinations) can be approximated as a finite weighted sum of exponentials [41], as given by $g_2=\sum _{i=0}^N (1/\theta _i) e^{a_i \tau }$. We solve Eq. (3) to obtain $I(t)$, for every $0 \leq i \leq n$, such that with mean $\mu =\tilde {\mu }/N$ and auto-correlation of the form $e^{a_i\tau }$. The resulting $I(t)'s$ were scaled by $1/\theta _i$ and added to get the final intensity of speckle pattern. Here we perform a non-linear least-squares fit to determine $a_i^{'}s$ for a fixed set of $\theta _i^{'}s$, where $a_i$ is the co-efficient of the fit, $\theta _i$ is the scaled factor and $n$ is the number of terms ($n=4$ is selected and $\theta _i$ were predetermined for these simulations). The terms $\tilde {\mu }$ and $\mu$ denote the mean intensity used in individual simulations (for every term $'n'$) and desired mean intensity. It is to be noted here that the individual simulations gives a random variables of the negative exponential distribution, while the resultant summation of random variable gives rise to gamma distribution [33,42] along with the desired complex autocorrelation.

While generating speckle intensity with the sum of such exponentials, the SDE in Eq. (3) can be modelled as, $dI(t)= -\alpha (I(t)- \alpha.p)dt +\sqrt {(2\alpha \mu I(t))} dW(t)$, which generates the speckles $I(t)$, with a PDF of the form $\frac {(1/\mu )^p I^{p-1} e^{-I/\mu }}{(J(p))}$; $J(p)=(p-1)\!$ and auto-correlation of the form $e^{-\alpha t}$. Note that the above PDF is a gamma distribution and a special case with $p=1$ leads to an exponential distribution. Here $p=1,2{\ldots }n,$ with $p$ being a positive integer. The Feller’s condition (i.e. $\frac {2\alpha \mu }{\sigma ^2} \geq 1$) for the above equation can be written as $\frac {2\alpha \mu p}{(\sqrt {2\alpha \mu })^2} = p$. As $p$ is a positive integer and $p\geq 1$, Feller’s condition [31] is always satisfied, thus generating positive speckle intensity for complex models.

2.3.3 Deep tissue blood flow

For deep tissue blood flow, we use a semi-infinite medium solution of the correlation diffusion equation (CDE) [2,37]. The above mentioned method for the generation of complex autocorrelation is used for the generation of the speckles. For the simulation, $\mu _a=0.1 cm^{-1}$, $\mu _s^{'}=12 cm^{-1}$ and $D_B=10^{-8} cm^{2}/s$ was used. In many cases, it may be sufficient to use a binomial approximation to the complex autocorrelation instead of adopting the approximation using weighted sum of simple exponential functions. The complex autocorrelation based on the infinite solution of CDE may be approximated (based on Binomial Approximation) to $g_1=e^{-\alpha _D \tau }$, where $\alpha _{D} =\frac {6 \mu _{s}^{'^2} k_{0}^2 D_{B}r}{2\sqrt {3\mu _{a}\mu _{s}^{'} }}$. Here, $\mu _a,\mu _s^{'},k_0$ and $D_B$ denote the absorption coefficient, reduced scattering coefficient, wave number and particle diffusion coefficient respectively. Although the results based on binomial approximation is not shown in this paper, it can be employed when the error in approximation is within the required accuracy.

2.3.4 In vivo animal experiment

To validate the SDE based algorithm, we have performed an in-vivo animal experiment in mice (Balp/C, 36 g, M, 4 months). This work was undertaken with the approval of the Institute Animal Ethical Committee (IAEC) at the APT research centre, Pune, India. The animal was anaesthetized using Ketamine and Xylazine and the skull of the mice was removed. For the uniform illumination, the laser beam (785 nm) was illuminated via a diffuser and for obtaining a point source (diameter $\leq 0.3 mm$), a focusing lens was used, as shown in Fig S1(a-b) in supplemental document. The detector used was a CCD camera (Basler acA 640-120 um) and the scattered intensity was measured at multiple exposure times. Using our recent multi-speckle DCS algorithm [16,43,44], the auto-correlation was measured and fitted to obtain the blood flow at an SD separation of 0.4 cm. The $\mu _a$ and $\mu _s^{'}$ were fixed as $0.1 cm^{-1}$ and $12 cm^{-1}$, respectively [2,45].

2.3.5 In vivo human experiment

A DCS system comprising a 785 nm laser (CW Thorlabs, L785P90) coupled to fibre that delivers light to the human forehead. At an SD separation of 2 cm, a single-mode fibre coupled to a single photon avalanche diode (SPAD, SPCM-AQRH-14 FC, Excelitas technologies) and a hardware autocorrelator (correlator.com) measures the intensity autocorrelation. The autocorrelation measured was used as the input to the SDE algorithm to generate the speckle intensity. The $\mu _a$ and $\mu _s^{'}$ were fixed as $0.12 cm^{-1}$ and $12 cm^{-1}$ respectively and the $D_B$ value obtained by fitting to the semi-infinite solution of CDE was $7.2 \times 10^{-9} cm^2/s$.

2.4 Noise model for speckle intensity pattern

Usually, for the deep tissue blood flow imaging (as in DCS), the auto-correlation is directly simulated either using Monte-Carlo simulation [46] or by solving CDE [37]. Instead, the proposed method simulates the speckle intensity (I(t)) using SDE. Hence, the appropriate noise model used in this work is the noise associated with the detection of intensity. For this, we use a camera-based noise model as follows: The noise model takes into account of quantum efficiency - QE, full well capacity - FWC, number of bits - n, pixel size - ps, dark current - $I_d$, exposure time - T, gain - g for a given sensor and irradiance - $I_r$. The model accounts for photon shot noise, dark noise, dark current shot noise and analogue to digital conversion (ADC) approximation based on Ref. [47]. The pseudocode for the same is given in Algorithm S1 of Supplemental document .

For the camera as a detector, the following specifications were used based on our recent work Ref: [16,45], Photometrics Prime BSI, with $QE-72\%$, pixel size $ps=5.6 \mu m$, $I_d$=68.6 $e^-$, FWC=45000$e^-$ and n=16). When an APD was used as a detector, the following parameters were used, QE=20$\%$, $Id=10^{-8} e^-$, gain=50. For APD-based detection, we did not integrate the speckle intensity for a given exposure time, and there was no quantization or thresholding based on FWC. For both camera and APD, note that we have predominantly used effects of quantum efficiency, dark noise, shot noise and gain of the detector, while we have neglected other types of noise such as readout noise, trans-impedance noise, fixed pattern noise, and source follower noise.

3. Results

3.1 Simulation of laser speckles for surface blood flow imaging

The simulated speckle intensity with the SDE algorithm for superficial blood flow are shown in Fig. 1. The sample is assumed to be uniformly illuminated by a laser source with a flow which gives rise to dynamic speckles with spatially independent simple exponential autocorrelation [1] given as $e^{-t/\tau _c}$, where $t$ is the correlation delay time and $\tau _c$ is the characteristics or decorrelation decay time. A sample discretised to $100 \times 100$ grid (representing pixels) with a high flow of $\tau _c=10^{-5} s$ is simulated in the centre channel and a $\tau _c=10^{-3} s$ is simulated on the sides, as shown in Fig. 1(a). We apply the SDE algorithm to every pixel to generate the speckle intensity both as a function of spatial variable (pixels) and time, as shown in Figs. 1(b) and 1(c), respectively. To confirm that the generated speckles have predetermined PDF and autocorrelation, we plot the histogram and intensity autocorrelation in Figs. 1(d) and (e), respectively. It is evident from the plots that the histogram fits a negative exponential PDF with an appropriate mean, and the autocorrelation fits a simple exponential decay with appropriate $\tau _c$ as desired. Here, the normalized intensity auto-correlation $g_2$ was computed based on the Siegert’s relation $g_2=1+\beta |g_1^2|$, where $g_1$ is the normalized field autocorrelation function and $\beta =1$ is the instrumentation constant. To confirm that every pixel obeys the desired autocorrelation, we fit the intensity autocorrelation $g_2$ for $\tau _c$ as shown in Fig. 1(g) and its inset figure, which agrees with the actual $\tau _c$ as shown in Fig. 1(f). The speckle intensity integrated to an exposure time of 1 ms are shown in Visualization 1 of the supplementary section as a video. To check the repeatability, we simulated speckles for 50 times, for a predetermined flow and mean intensity, which resulted in a mean percentage error of $2.8\%$ and $3\%$ for $\tau _c$ and mean intensity respectively. In addition to auto-correlation, we also validate the intensity of speckles in terms of speckle contrast. Figure 1(h) and Fig. 1(i) shows the speckle contrast produced spatially [1] (using a $3 \times 3$ window) and temporally [12] (over 100 frames), respectively, for the data generated at $100 \mu s$ exposure time, both indicating the speckle contrast has lower values when the flow is high at the centre.

 figure: Fig. 1.

Fig. 1. Simulation of speckle intensity using SDE algorithm for surface blood flow imaging. (a) Schematic diagram of the LSCI system simulated by the SDE algorithm, wherein a laser source (LD) is uniformly illuminated over the sample, and the camera captures the scattered intensity. The sample has a high flow region at the centre bounded by low flow regimes. (b) The speckle intensity simulated using the proposed method are shown (Refer to the supplementary video Visualization 1 for the intensity speckle fluctuations as a function of time and space), wherein the sample is discretized into $100 \times 100$ pixels. The model used for the desired auto-correlation is $e^{-t/\tau _c}$, wherein the $\tau _c$ for the low and high flow regions are $10^{-5}$s and $10^{-3}$s respectively. (c) The simulated speckle intensity for a given pixel is plotted as a function of time, t, where the inset picture shows the zoomed version of the same. (d) The histogram of the speckle generated for a given mean $\mu$ by the proposed method produces a PDF with negative exponential distribution as desired. (e) The temporal autocorrelation of simulated speckles and the desired speckles. The desired flow is represented in terms of $\tau _c$ is shown in (f), whereas (g) shows the fitted $\tau _c$ from the simulated speckles. (h) and (i) respectively show the spatial and temporal speckle contrast obtained from the speckle simulated at an exposure time of 1 ms. The spatial speckle contrast was obtained using a window of $3 \times 3$ and temporal speckle contrast was obtained using 100 images.

Download Full Size | PDF

It is well known that the sum of random variables with negative exponential distributions gives rise to a resultant random variable with gamma distribution [33]. The same concept applies to speckles where we sum up several speckles to get a resultant speckle with gamma distribution which is also associated with a decrement of instrumentation factor $\beta$ [26] to be less than one. In such a scenario, the autocorrelation of the resultant speckles will be a weighted sum of individual autocorrelations (of component speckles) [38]. Using the SDE algorithm, we simulate several speckle intensity and perform a summation of speckles, as shown in Fig. 2. The decrement in the instrumentation factor $\beta$ also summarizes the importance of the appropriate speckle to pixel mapping in the laser speckle based blood flow measurement systems, and our results of decrement in $\beta$ when the speckles from multiple pixels are summed up is consistent with earlier literature [26].

 figure: Fig. 2.

Fig. 2. Generating speckles with more complex auto-correlation. This figure illustrates the addition of the speckles generated by the proposed SDE algorithm (each following negative exponential PDF and a simple exponential auto-correlation) for four independent pixels. The resultant speckle follows gamma distribution and the desired auto-correlation.

Download Full Size | PDF

3.1.1 Simulation of different flow models for various tissue types in LSCI

Different flow models having more complex field autocorrelation (in contrast to simple exponential relations) that may arise in different parts of the tissues, as explained in Ref [6], are simulated by the SDE method. Using the above-mentioned principles based on the addition of speckles, as depicted in Fig. 2, the complex auto-correlation models (as given in Eq. (5) in the Methods section) were generated for different cases, such as with and without static layers [29], with different $n$ and $d$ values based on the tissue type. The autocorrelation of the speckle intensity generated using SDE is compared against the desired auto-correlation models as shown in Fig. 3 with a reasonably good fit, thus indicating the capability of the proposed method to model complex blood flow models.

 figure: Fig. 3.

Fig. 3. Simulation of complex auto-correlation for different tissue types. For three representative tissue types such as large vessels, small vessels, and parenchyma’s, the complex autocorrelation models based on different types of the dynamic scattering regimes $d$, amount of static scatterers $\rho$ and $\tau _c$ were generated. These were used for the generation of speckles using SDE and their auto-correlation was compared with the desired auto-correlation.

Download Full Size | PDF

3.1.2 In-vivo experimental validation of laser speckle contrast imaging

For in-vivo validation of the proposed algorithm, we measured the surface blood flow in the mice brain using our recently proposed Multi-speckle imaging system [16,43] as shown in Fig. 4(a). The system uses a low-frame rate camera to measure the speckle intensity at multiple exposure times and compute the intensity autocorrelation, which is related to blood flow. The schematic of the experimental setup is shown in Fig S1 (a) of the supplemental section. For the simulation of speckles, the intensity image of the region of interest in the brain at 1.1 ms exposure time was used as the mean intensity at every pixel. The characteristic decay time $\tau _c$ determined experimentally is shown in Fig. 4(b). For two representative pixels, the multi exposure speckle contrast data obtained by the experiment is compared with the simulated speckle contrast data. The corresponding $\tau _c$ values are reported in the legend of inset of Fig. 4(b) and are comparable to one another. The SDE based intensity data at $1.1 ms$ exposure time is shown in Fig. 4(c). Visualization 2 (Supplemental section) shows the simulated speckle intensity fluctuations in mice brain. The temporal speckle contrast obtained from the simulated SDE-based speckle intensity is shown for a given exposure time of $T=1.1 ms$ is compared with the experimental data.

 figure: Fig. 4.

Fig. 4. Experimental validation of SDE-based simulation of laser speckles. (a) Illustration of the LSCI experimental system used for imaging mice brain. (b) The experimental and the simulated SDE-based speckle intensity were fitted for $e^{-t/\tau _c}$ model to obtain the decorrelation time $\tau _c$. (c) The simulated speckle intensity is integrated upto 1 ms exposure time of the camera. (d) The Experimental and simulated speckle contrast from the simulated speckle intensity at exposure time T=1.1 ms.

Download Full Size | PDF

3.2 Simulation of laser speckles for deep tissue blood flow imaging

The schematic of the DCS system with a pointed source used for simulation is shown in Fig. 5(a). A laser source is pointed to the centre of the object, and the light diffuses symmetrically in all directions. The sample is assumed to be homogenous, and the optical properties used for the simulation are $\mu _a=0.1 cm^{-1}$, $\mu _s^{'}=0.12 cm^{-1}$ and the dynamic property (particle diffusion coefficient) used is $D_B=10^{-8} cm^2/s$. An infinite solution of Diffusion equation [2,48] and semi-infinite CDE [2,37] are used to generate the intensity and autocorrelation, respectively, for the simulated speckles. The SDE-based algorithm was used to simulate the speckle intensity and the integrated intensity $(\phi )$ at an exposure time of T=1ms as shown in Fig. 5(c), which is comparable to the desired intensity as shown in Fig. 5(b). The log of the intensity is represented in Fig. 5(b, c, e) for better clarity. The intensity profile as a function of SD separation $(log(\phi.r))$ is shown in Fig. 5(d) along with a linear fit. It is well known that the slope of the linear fit should approximately follow $k=\sqrt {3\mu _a \mu _s^{'}} =2.08$, which is the desired value, whereas the slope obtained from the SDE-based intensity speckle was 2.1. The speckle intensity at an exposure time of $10 \mu s$ is shown in Fig. 5(e) (Refer to the supplementary video Visualization 3 for the intensity of speckle fluctuations as a function of time and space). To validate the proposed method, we have compared the histogram (PDF) and the autocorrelation of the speckles obtained by the simulation against their desired values. Figure 5(f) compares the SDE-based intensity autocorrelation at two representative SD separations (r=1cm and r=3 cm) against the desired auto-correlation, which are in reasonable comparison with one another. Figure 5(g) shows the histogram of the speckle intensity obtained at a given representative SD separation which fits a gamma distribution with the desired scale parameter (Desired $\mu =2.26 \times 10^{-15}$ and the fitted $\mu =1.85 \times 10^{-15}$). Figures 5(h) and 5(i) shows the intensity auto-correlations obtained at two delay times $\tau =10 \mu s$ and $\tau =100 \mu s$ respectively. In addition, we also plot the computed speckle contrast as a function of multiple exposure times for two different SD separations ($r=1 cm$ and $3 cm$) in Fig. 5(j). These results show that the SDE algorithm can be used for the generation of both uniform as well as diffusing speckle intensity with applications in conventional laser speckle imaging and deep-tissue imaging, respectively.

 figure: Fig. 5.

Fig. 5. Simulation of laser speckles for deep tissue blood flow. Schematic diagram of the DCS system simulated by the SDE algorithm, wherein a laser source (LD) is focused on the sample, and the detector captures the scattered intensity. The sample had the following optical and dynamic properties: $\mu _a=0.1 cm^{-1},\mu _s^{'}=0.12 cm^{-1} D_B=10^{-8} cm^2/s$. (b) The desired intensity has a diffusing profile from the centre (of the source) towards the boundary. (c) The simulated speckle intensity using the proposed method are shown at an exposure time of 1 ms. Both (b) and (c) show the log of the speckle intensity for better visualization. (d) The profile of the simulated speckles (in terms of $log(\phi.r)$) is fitted to a linear plot with the slope comparable to the desired slope indicating the match of optical properties. (e) The simulated speckle intensity integrated up to 10 $\mu s$ exposure time. (Refer to the supplementary video Visualization 3 for the intensity speckle fluctuations as a function of time and space). (f) The normalized intensity auto-correlation $g_2$ at two representative SD pairs, r =1 cm and 3 cm, along with the desired auto-correlation obtained using the semi-infinite solution of CDE. (g) The histogram of the simulated intensity speckle fits to a gamma distribution with the desired $\mu$. (h) and (i) shows the auto-correlation $g_2$ at delay times $\tau =10 \mu s$ and $100 \mu s$ as a function of spatial coordinates. (j) Speckle contrast is plotted for multiple exposure times at two representative SD pairs (r=1 cm and 3 cm).

Download Full Size | PDF

3.2.1 In-vivo experimental validation of diffusing speckles generation

The capability of the SDE algorithm to generate speckle intensity comparable to $in-vivo$ experimental studies is shown in Fig. 6. The results show the in-vivo validation performed in both mice and adult human brain for the speckles generation using the SDE algorithm. The experimental method is detailed in the methods section. Figure 6(a) shows the comparison between the experimental intensity autocorrelation at an SD separation of 0.4 cm in the brain of a mouse in alive and sacrificed condition against the SDE-based intensity autocorrelation. The $g_2$ curves obtained experimentally were fitted to the semi-infinite solution of CDE, which resulted in $D_B$ of $1.57 \times 10^{-9} cm^2/s$ and $1.49 \times 10^{-11} cm^2/s$ in alive and sacrificed condition. The auto-correlation of speckle intensity generated by SDE resulted in $D_B$ of $1.45 \times 10^{-9} cm^2/s$ and $1.53 \times 10^{-11} cm^2/s$. Figure 6(b) shows the experimental autocorrelation in an adult human head at an SD separation of 2 cm compared to the SDE-based intensity speckle simulation. These results indicate that SDE-based algorithms can be used to obtain speckle intensity for deep tissue blood flow comparable with the experimental results based on $in-vivo$ animal and human studies.

 figure: Fig. 6.

Fig. 6. In-vivo animal and human experimental results for validating the proposed algorithm. (a) The autocorrelation $g_2$ of the experimental and simulated diffuse speckle intensity measured at an SD separation of 0.4cm in the brain of a mouse in alive and sacrificed condition. (b) The DCS measurement was performed in an adult human brain at an SD of 2 cm, and the normalized field auto-correlation was extracted based on Siegert’s relation. Intensity autocorrelation obtained from the simulated speckle is comparable to the experimental results.

Download Full Size | PDF

3.3 Modelling of speckle intensity with realistic noise models

The speckle intensity obtained using SDE were further added with detector noise for more realistic analysis. For this, we assume that our detectors were (a) a camera and (b) an APD [2]. If the detector is a camera, we use the camera noise model for the speckle intensity spatially. For the APD, we apply the noise to temporal intensity fluctuations at a given pixel. Figure 7(a) shows the speckle intensity with a camera noise model with dark and shot noise effects. The speckle intensity deviated from the linear path when dark noise was added. Figure 7(b) shows the effect of the noise on an APD, wherein the noise was added to the speckle intensity for a given SD pair. It can be seen that SDE-based autocorrelation without noise follows the desired autocorrelation, while speckle intensity with added noise has the autocorrelation deviating slightly from the desired value. It is also to be noted that the effect of noise is evident at larger values of delay time $\tau$ than that of earlier times as the number of samples available is large at earlier delay times which may results in averaging of noise. The effect of the noises added to speckle intensity is more realistic than that of the noise models used earlier in DCS, as the DCS usually follows the noise models for auto-correlation, while the noises should be applied to the fundamental quantity which is speckle intensity. Figure 7(c) shows the camera noise effects on normalized speckle contrast $(\kappa _N)$. The speckle contrast fits the initial exposure times; however, as the exposure time increases, due to inadequate dynamic range or full-well capacity of the camera, it saturates, which is evident from the plot.

 figure: Fig. 7.

Fig. 7. Effect of noise added to speckle intensity for realistic modelling of speckle intensity. (a) Intensity profile $log(\phi.r)$) is plotted as a function of $r$ at an exposure time of $1\mu s$. When no noise is added, the SDE generates an intensity profile that follows a linear fit. When the camera noise was added, the SNR decreases as the SD separation increases and the signal deviates from the linear fit. (b) At a given representative SD separation, speckle intensity generated by SDE follows the desired autocorrelation when no noise is added. The detector noise is applied (based on APD), resulting in the deviation of the autocorrelation from the desired value for larger values of delay times. (c) The simulated speckle contrast at various exposure times $T$ is plotted. At larger time $T$, due to saturation of the camera, the speckle contrast deviates from the desired value. (d) The mean and the variance of the multi-exposure speckle contrast calculated by varying the number of speckles $N$. As the number of speckles used to calculate speckle contrast increases, the standard deviation decreases and mean speckle contrast approaches the analytical solution.

Download Full Size | PDF

Although we have used a $3 \times 3$ window for obtaining spatial speckle contrast as shown in Fig. 1(h); as reported in Ref [49], spatial window of different sizes may used for optimal computation of speckle contrast. In case of calculation of temporal speckle contrast as shown in Fig. 1(i), it is well known that, more the number of pixels/ frames used, better would be the SNR of the speckle contrast, as shown in Ref [50] in the noise model for speckle contrast. Figure 7(d) shows the mean and the variance of the multi-exposure speckle contrast calculated by varying the number of speckles. As the number of speckles used to calculate speckle contrast increases, the standard deviation decreases and mean speckle contrast approaches the analytical solution.

4. Discussion

This paper proposes a comprehensive tool for laser intensity speckle simulation that uses a novel stochastic differential equations-based algorithm. The proposed algorithm generates speckle intensity that follow a desired PDF and autocorrelation function for both surface and deep tissue blood flow. The speckle intensity produced by SDE based algorithm were validated by extensive simulation studies and in-vivo experiments measuring speckles from mice brain and an adult human head.

For the surface blood flow, we simulate the speckle intensity with a simple field autocorrelation given by $g_1=exp(-t/\tau _c)$ for different values of $\tau _c$ as given in Fig. 1. We validate the SDE algorithm by comparing the two main characteristics of the speckle intensity [33]: the PDF of the speckles and autocorrelation, against their respective desired values as shown in Fig. 1(d,e). In addition, we generate speckle contrast [1] (Fig. 1(h)), which is one of the widely used quantifiers in measuring surface blood flow. We also validate the LSCI system for different tissue types, such as large vessels, small vessels and parenchyma types for different values of static scatterers and different orders of the diffusing regime based on different flow models. This shows that SDE could be a potential tool for the simulation of speckle intensity for different flow models, as it is established in Ref [6] that different tissue types follow different flow models in measuring surface blood flow. An in-vivo validation of the proposed SDE algorithm in mice brain is shown in Fig. 4 by comparing simulated speckles against the in-vivo experimental speckle measurements. Thus, the SDE algorithms could be used as a robust tool for generating speckle intensity for LSCI with different flow models.

In order to simulate the deep tissue blood flow speckles, we combine the SDE with the solution of the photon diffusion equation along with appropriate PDF and autocorrelation. The infinite solution of the photon diffusion equation is used for the generation of the mean intensity (or initial intensity $\mu$ in Eq. (3)). Similarly, the semi-infinite solution of CDE was used to select the normalized field autocorrelation function [2], which is not a simple autocorrelation function of exponential type. In order to accommodate complicated (not simple exponential) autocorrelation, we adopt the method of approximating such functions using weighted series of simple exponentials. The speckles simulated for each such simple exponentials are appropriately added to get the speckles with desired PDF and autocorrelation. It can be deduced that such speckles obey gamma distribution, which is the case for the sum of several spatially averaged speckles due to the finite size of the detector [50,51]. The simulation results as shown in Fig. 5 and the in-vivo results on mice and healthy human adults as shown in Fig. 6 validate the proposed algorithm. Additionally, the multi-exposure speckle contrast, a derived quantity used by several modalities was also generated from the SDE based speckle intensity and its validation is shown in Fig. 5(k).

In order to simulate speckles in a more realistic manner, we have also incorporated the noise models associated with commonly employed detectors for speckle measurement. The noise model considers critical detector parameters such as quantum efficiency, dark noise, pixel size / active detection area, gain, full well capacity, dynamic range, and quantization noise. The proposed noise models are directly associated with the speckle intensity rather than the existing noise models that have been developed for either the autocorrelation or the speckle contrast, which are derived quantities from the speckle as in LSCI/DCS [50,52]. We have not included the effects laser coherence length, line width of the laser and numerical aperture of the camera lens in the model. It is shown in Ref [53] that, volume holographic grating stabilized laser diode and other laser sources whose spectral width is less than 1 nm improves the speckle contrast. In addition, a complete theoretical study on speckle contrast and laser coherence length and imaging characteristics of the optical system is reported in Ref [54]. These models could be incorporated to model the effect of coherence length and line width in order to access its effect on the speckle intensity produced by the proposed method.

As mentioned in Section 2.3, we have used $t_{min}=10^{-6} s,t_{max}=1s$ and $N=10^7$ samples to generate the speckle intensity. In this time discretization, within $5\%$ error in blood flow obtained using intensity autocorrelation, we can simulate the flow or $\tau _c$ in the range of $5\mu s$ to $40 ms$. In principle, the SDE algorithm could simulate speckle intensity without any constraints on blood flow. However, considering the practical computational constraints such as the numerical instability of the algorithm and the requirement of a sufficient number of samples for the calculation of intensity autocorrelation, the above range is found to be optimal. However, we can increase the velocity range by varying the time range of the simulations. Additionally, to check the SDE algorithm’s repeatability, we produced speckle intensity for a predetermined flow and mean intensity over 50 trails. We found that the mean error in the flow and mean intensity produced is around $2.8\%$ and $3\%$, respectively, which shows that the algorithm is repeatable. The noise characteristics of the camera, such as, sensor full well capacity, quantum efficiency, and signal-to-noise ratio, are added later after the simulation of speckle intensity. The above-mentioned parameters determines the usable range of flow thus mimicking the realistic case. For instance, for the given camera (Basler acA-640-120 um), based on its quantum efficiency, full-well capacity and noise characteristics, $\tau _c$ should be within the range of $0.1$ to $10 ms$, beyond which the intensity saturates due to the finite dynamic range of the camera. It is crucial to use the unbiased auto-correlation as given in the Methods section for auto-correlation calculation. In case of biased/scaled auto-correlation models are used, it could lead to bias in the larger values of $\tau$, as shown in Fig S3 in supplemental document. The bias is caused by the insufficient number of samples for the auto-correlation at larger values of $\tau$, whose effects are removed when an unbiased auto-correlation is used. Additionally, it is evident that we have generated temporal intensity speckles at every pixels and not spatially across multiple pixels. Under the assumption of ergodicity, the temporal speckle contrast and spatial speckle contrast are equivalent. However, further mathematical analysis, as given in Ref [42,55], would be required to ensure the same.

Although there are a few methods to produce speckle intensity, as given in Ref [26], most do not provide a means to apriori select complex (non-exponential) autocorrelation. Even those which could incorporate autocorrelation statistics in blood flow are usually a simple exponential function as in Ref [28,30,56]. We have introduced the complex autocorrelation models for both surface as well as deep tissue blood flow, by approximation using finite weighted exponential series. This allows us to model speckles from any tissues with multiple layers of optical and dynamical properties. This is essential, as the inappropriate forward modelling leads to significantly large errors in the reconstruction of blood flow [57,58] as encountered in modelling the flow in the human head and animal brain with multiple layers [59].

The proposed algorithm could simulate speckle intensity of any desired auto-correlation and speckle statistics by appropriately selecting the parameters $a$ and $b$ of Eq. (3). Although the proposed method in this paper, shows the applicability of producing speckle intensity for laser speckle contrast imaging, it could also be easily extended for other speckle statistics, which are not exponential or gamma distribution [60]. Moreover, although we have shown the commonly used reflection geometry for validation, simulations can also be implemented for transmission geometry which gives better signal-to-background ratio [61]. Additionally, in near future, the proposed laser simulation tool could be potentially used to compare multiple laser speckle imaging modalities and could optimize several parameters such as exposure time, number of frames/ speckles, optode positions, which are key factors for both LSCI and DCS based blood flow measurement systems.

The SDE-based speckle simulation for turbid medium like tissue is a new approach based on stochastic process whereas the most of the existing methods use statistical tools, ODE’s or Fourier-based approaches [2628]. We have validated the proposed method using extensive simulation and in-vivo studies, thus proving it to be a robust, comprehensive tool to simulate speckle intensity for both surface and deep tissue flow imaging applications. We believe this approach would eventually help adopt stochastic differential equations to model light propagation in a dynamic turbid medium like tissue.

Funding

Science and Engineering Research Board (Early Career Research Award -2016); Department of Biotechnology, Ministry of Science and Technology, India (Ramalingaswamy Fellowship-2016); Indian Institute of Technology Bombay (Seed Grant, Wadhwani Research Center for Bioengineering (WRCB)).

Disclosures

The authors have no conflict of interest.

Data availability

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

Supplemental document

See Supplement 1 for supporting content.

References

1. D. A. Boas and A. K. Dunn, “Laser speckle contrast imaging in biomedical optics,” J. Biomed. Opt. 15(1), 011109 (2010). [CrossRef]  

2. T. Durduran, R. Choe, W. Baker, and A. G. Yodh, “Diffuse optics for tissue monitoring and tomography,” Rep. Prog. Phys. 73(7), 076701 (2010). [CrossRef]  

3. J. Briers and A. Fercher, “Retinal blood-flow visualization by means of laser speckle photography,” Investigative ophthalmology & visual science 22(2), 255–259 (1982).

4. A. B. Parthasarathy, W. J. Tom, A. Gopal, X. Zhang, and A. K. Dunn, “Robust flow measurement with multi-exposure speckle imaging,” Opt. Express 16(3), 1975–1989 (2008). [CrossRef]  

5. W. Liu, R. Qian, S. Xu, P. Chandra Konda, J. Jönsson, M. Harfouche, D. Borycki, C. Cooke, E. Berrocal, Q. Dai, H. Wang, and R. Horstmeyer, “Fast and sensitive diffuse correlation spectroscopy with highly parallelized single photon detection,” APL Photonics 6(2), 026106 (2021). [CrossRef]  

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

7. J. D. Johansson, D. Portaluppi, M. Buttafava, and F. Villa, “A multipixel diffuse correlation spectroscopy system based on a single photon avalanche diode array,” J. Biophotonics 12(11), e201900091 (2019). [CrossRef]  

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

9. T. Durduran and A. G. Yodh, “Diffuse correlation spectroscopy for non-invasive, micro-vascular cerebral blood flow measurement,” NeuroImage 85, 51–63 (2014). [CrossRef]  

10. J. P. Culver, T. Durduran, D. Furuya, C. Cheung, J. H. Greenberg, and A. Yodh, “Diffuse optical tomography of cerebral blood flow, oxygenation, and metabolism in rat during focal ischemia,” J. Cereb. Blood Flow Metab. 23(8), 911–924 (2003). [CrossRef]  

11. R. Bi, J. Dong, and K. Lee, “Deep tissue flowmetry based on diffuse speckle contrast analysis,” Opt. Lett. 38(9), 1401–1403 (2013). [CrossRef]  

12. C. P. Valdes, H. M. Varma, A. K. Kristoffersen, T. Dragojevic, J. P. Culver, and T. Durduran, “Speckle contrast optical spectroscopy, a non-invasive, diffuse optical method for measuring microvascular blood flow in tissue,” Biomed. Opt. Express 5(8), 2769–2784 (2014). [CrossRef]  

13. H. M. Varma, C. P. Valdes, A. K. Kristoffersen, J. P. Culver, and T. Durduran, “Speckle contrast optical tomography: A new method for deep tissue three-dimensional tomography of blood flow,” Biomed. Opt. Express 5(4), 1275–1289 (2014). [CrossRef]  

14. W. Zhou, O. Kholiqov, S. P. Chong, and V. J. Srinivasan, “Highly parallel, interferometric diffusing wave spectroscopy for monitoring cerebral blood flow dynamics,” Optica 5(5), 518–527 (2018). [CrossRef]  

15. C. Huang, D. Irwin, Y. Lin, Y. Shang, L. He, W. Kong, J. Luo, and G. Yu, “Speckle contrast diffuse correlation tomography of complex turbid medium flow,” Med. Phys. 42(7), 4000–4006 (2015). [CrossRef]  

16. K. Murali and H. M. Varma, “Multi-speckle diffuse correlation spectroscopy to measure cerebral blood flow,” Biomed. Opt. Express 11(11), 6699–6709 (2020). [CrossRef]  

17. R. Paul, K. Murali, S. Chetia, and H. M. Varma, “High-density diffuse correlation tomography with enhanced depth localization and minimal surface artefacts,” Biomed. Opt. Express 13(11), 6081–6099 (2022). [CrossRef]  

18. W. Zhou, M. Zhao, O. Kholiqov, and V. J. Srinivasan, “Multi-exposure interferometric diffusing wave spectroscopy,” Opt. Lett. 46(18), 4498–4501 (2021). [CrossRef]  

19. J. Xu, A. K. Jahromi, J. Brake, J. E. Robinson, and C. Yang, “Interferometric speckle visibility spectroscopy (isvs) for human cerebral blood flow monitoring,” APL Photonics 5(12), 126102 (2020). [CrossRef]  

20. A. Biswas, A. Takshi, and A. B. Parthasarathy, “A method for improving the dynamic range of integrated diffuse speckle contrast spectroscopy,” in Optics and the Brain, (Optical Society of America, 2021), pp. BTh1B–4.

21. E. James and S. Powell, “Fourier domain diffuse correlation spectroscopy with heterodyne holographic detection,” Biomed. Opt. Express 11(11), 6755–6779 (2020). [CrossRef]  

22. M. Seong, Y. Oh, K. Lee, and J. G. Kim, “Blood flow estimation via numerical integration of temporal autocorrelation function in diffuse correlation spectroscopy,” Comput. Methods Programs Biomed. 222, 106933 (2022). [CrossRef]  

23. Z. Hajjarian, H. T. Nia, S. Ahn, A. J. Grodzinsky, R. K. Jain, and S. K. Nadkarni, “Laser speckle rheology for evaluating the viscoelastic properties of hydrogel scaffolds,” Sci. Rep. 6(1), 37949 (2016). [CrossRef]  

24. H. M. Varma, K. P. Mohanan, N. Hyvönen, A. K. Nandakumaran, and R. M. Vasu, “Ultrasound-modulated optical tomography: recovery of amplitude of vibration in the insonified region from boundary measurement of light correlation,” J. Opt. Soc. Am. A 28(11), 2322–2331 (2011). [CrossRef]  

25. R. S. Chandran, S. Sarkar, R. Kanhirodan, D. Roy, and R. M. Vasu, “Diffusing-wave spectroscopy in an inhomogeneous object: local viscoelastic spectra from ultrasound-assisted measurement of correlation decay arising from the ultrasound focal volume,” Phys. Rev. E 90(1), 012303 (2014). [CrossRef]  

26. S. J. Kirkpatrick, D. D. Duncan, and E. M. Wells-Gray, “Detrimental effects of speckle-pixel size matching in laser speckle contrast imaging,” Opt. Lett. 33(24), 2886–2888 (2008). [CrossRef]  

27. D. D. Duncan and S. J. Kirkpatrick, “The copula: a tool for simulating speckle dynamics,” J. Opt. Soc. Am. A 25(1), 231–237 (2008). [CrossRef]  

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

29. H. J. Rabal and R. A. Braga Jr, Dynamic Laser Speckle and Applications (CRC press, 2018).

30. E. James, S. Powell, and P. Munro, “Simulation of statistically accurate time-integrated dynamic speckle patterns in biomedical optics,” Opt. Lett. 46(17), 4390–4393 (2021). [CrossRef]  

31. C. A. Braumann, Introduction to Stochastic Differential Equations with Applications to Modelling in Biology and Finance (John Wiley & Sons, 2019).

32. C. Gardiner, “Stochastic methods: a handbook for the natural and social sciences 4th ed. (2009),”.

33. J. W. Goodman, Speckle Phenomena in Optics: Theory and Applications (Roberts and Company Publishers, 2007).

34. G. N. Mil’shtein, “Approximate integration of stochastic differential equations,” Teoriya veroyatnostei i ee primeneniya 19, 583–588 (1974).

35. G. Mil’shtejn, “Approximate integration of stochastic differential equations,” Theory Probab. & Its Appl. 19(3), 557–562 (1975). [CrossRef]  

36. Z. Yin and S. Gan, “An improved milstein method for stiff stochastic differential equations,” Adv. Differ. Eqs. 2015(1), 369 (2015). [CrossRef]  

37. D. A. Boas, L. Campbell, and A. G. Yodh, “Scattering and imaging with diffusing temporal field correlations,” Phys. Rev. Lett. 75(9), 1855–1858 (1995). [CrossRef]  

38. R. Zárate-Mi nano and F. Milano, “Construction of sde-based wind speed models with exponentially decaying autocorrelation,” Renewable Energy 94, 186–196 (2016). [CrossRef]  

39. D. J. Higham, “An algorithmic introduction to numerical simulation of stochastic differential equations,” SIAM Rev. 43(3), 525–546 (2001). [CrossRef]  

40. E. M. Buckley, A. B. Parthasarathy, P. E. Grant, A. G. Yodh, and M. A. Franceschini, “Diffuse correlation spectroscopy for measurement of cerebral blood flow: future prospects,” Neurophotonics 1(1), 011009 (2014). [CrossRef]  

41. B. M. Bibby, I. M. Skovgaard, and M. Sørensen, “Diffusion-type models with given marginal distribution and autocorrelation function,” Bernoulli 11(2), 191–220 (2005). [CrossRef]  

42. O. Ibe, Markov Processes for Stochastic Modeling (Newnes, 2013), Chapter 1, 2 and chapter 9.

43. K. Murali, A. Nandakumaran, T. Durduran, and H. M. Varma, “Recovery of the diffuse correlation spectroscopy data-type from speckle contrast measurements: towards low-cost, deep-tissue blood flow measurements,” Biomed. Opt. Express 10(10), 5395–5413 (2019). [CrossRef]  

44. K. Murali, A. Nandakumaran, and H. M. Varma, “On the equivalence of speckle contrast-based and diffuse correlation spectroscopy methods in measuring in vivo blood flow,” Opt. Lett. 45(14), 3993–3996 (2020). [CrossRef]  

45. R. Paul, K. Murali, S. Chetia, and H. M. Varma, “A simple algorithm for diffuse optical tomography without Jacobian inversion,” Biomed. Phys. Eng. Express 8(4), 045001 (2022). [CrossRef]  

46. Q. Fang and D. A. Boas, “Monte carlo simulation of photon migration in 3d turbid media accelerated by graphics processing units,” Opt. Express 17(22), 20178–20190 (2009). [CrossRef]  

47. M. Konnik and J. Welsh, “High-level numerical simulations of noise in ccd and cmos photosensors: review and tutorial,” arXiv, arXiv:1412.4031 (2014). [CrossRef]  

48. S. R. Arridge, M. Cope, and D. Delpy, “The theoretical basis for the determination of optical pathlengths in tissue: temporal and frequency analysis,” Phys. Med. Biol. 37(7), 1531–1560 (1992). [CrossRef]  

49. D. D. Duncan, S. J. Kirkpatrick, and R. K. Wang, “Statistics of local speckle contrast,” J. Opt. Soc. Am. A 25(1), 9–15 (2008). [CrossRef]  

50. S. Yuan, “Sensitivity, noise and quantitative model of laser speckle contrast imaging,” Ph.D. thesis, Tufts University (2008).

51. J. C. Dainty, Laser Speckle and Related Phenomena, vol. 9 (Springer science & business Media, 2013).

52. C. Zhou, G. Yu, D. Furuya, J. H. Greenberg, A. G. Yodh, and T. Durduran, “Diffuse optical correlation tomography of cerebral blood flow during cortical spreading depression in rat brain,” Opt. Express 14(3), 1125–1144 (2006). [CrossRef]  

53. D. D. Postnov, X. Cheng, S. E. Erdener, and D. A. Boas, “Choosing a laser for laser speckle contrast imaging,” Sci. Rep. 9(1), 2542 (2019). [CrossRef]  

54. H. Fujii and T. Asakura, “A contrast variation of image speckle intensity under illumination of partially coherent light,” Opt. Commun. 12(1), 32–38 (1974). [CrossRef]  

55. E. Parzen, “Conditions that a stochastic process be ergodic,” Ann. Math. Statist. 29(1), 299–301 (1958). [CrossRef]  

56. X. Cheng, H. Chen, E. J. Sie, F. Marsili, and D. A. Boas, “Development of a monte carlo-wave model to simulate time domain diffuse correlation spectroscopy measurements from first principles,” J. Biomed. Opt. 27(08), 083009 (2022). [CrossRef]  

57. L. Gagnon, M. Desjardins, J. Jehanne-Lacasse, L. Bherer, and F. Lesage, “Investigation of diffuse correlation spectroscopy in multi-layered media including the human head,” Opt. Express 16(20), 15514–15530 (2008). [CrossRef]  

58. H. Zhao, E. Sathialingam, and E. M. Buckley, “Accuracy of diffuse correlation spectroscopy measurements of cerebral blood flow when using a three-layer analytical model,” Biomed. Opt. Express 12(11), 7149–7161 (2021). [CrossRef]  

59. M. M. Wu, K. Perdue, S.-T. Chan, K. A. Stephens, B. Deng, M. A. Franceschini, and S. A. Carp, “Complete head cerebral sensitivity mapping for diffuse correlation spectroscopy using subject-specific magnetic resonance imaging models,” Biomed. Opt. Express 13(3), 1131–1151 (2022). [CrossRef]  

60. N. Bender, H. Yılmaz, Y. Bromberg, and H. Cao, “Customizing speckle intensity statistics,” Optica 5(5), 595–600 (2018). [CrossRef]  

61. D.-Y. Li, Q. Xia, T.-T. Yu, J.-T. Zhu, and D. Zhu, “Transmissive-detected laser speckle contrast imaging for blood flow monitoring in thick tissue: from monte carlo simulation to experimental demonstration,” Light: Sci. Appl. 10(1), 241 (2021). [CrossRef]  

Supplementary Material (4)

NameDescription
Supplement 1       Supplemental Document
Visualization 1       1: Simulation of intensity speckles using SDE algorithm for surface blood flow imaging
Visualization 2       2: Simulation of intensity speckles for in-vivo validation using SDE algorithm in mice brain (Near Bregma)
Visualization 3       3: Simulation of intensity speckles using SDE algorithm for deep tissue blood flow imaging

Data availability

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

Cited By

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

Alert me when this article is cited.


Figures (7)

Fig. 1.
Fig. 1. Simulation of speckle intensity using SDE algorithm for surface blood flow imaging. (a) Schematic diagram of the LSCI system simulated by the SDE algorithm, wherein a laser source (LD) is uniformly illuminated over the sample, and the camera captures the scattered intensity. The sample has a high flow region at the centre bounded by low flow regimes. (b) The speckle intensity simulated using the proposed method are shown (Refer to the supplementary video Visualization 1 for the intensity speckle fluctuations as a function of time and space), wherein the sample is discretized into $100 \times 100$ pixels. The model used for the desired auto-correlation is $e^{-t/\tau _c}$, wherein the $\tau _c$ for the low and high flow regions are $10^{-5}$s and $10^{-3}$s respectively. (c) The simulated speckle intensity for a given pixel is plotted as a function of time, t, where the inset picture shows the zoomed version of the same. (d) The histogram of the speckle generated for a given mean $\mu$ by the proposed method produces a PDF with negative exponential distribution as desired. (e) The temporal autocorrelation of simulated speckles and the desired speckles. The desired flow is represented in terms of $\tau _c$ is shown in (f), whereas (g) shows the fitted $\tau _c$ from the simulated speckles. (h) and (i) respectively show the spatial and temporal speckle contrast obtained from the speckle simulated at an exposure time of 1 ms. The spatial speckle contrast was obtained using a window of $3 \times 3$ and temporal speckle contrast was obtained using 100 images.
Fig. 2.
Fig. 2. Generating speckles with more complex auto-correlation. This figure illustrates the addition of the speckles generated by the proposed SDE algorithm (each following negative exponential PDF and a simple exponential auto-correlation) for four independent pixels. The resultant speckle follows gamma distribution and the desired auto-correlation.
Fig. 3.
Fig. 3. Simulation of complex auto-correlation for different tissue types. For three representative tissue types such as large vessels, small vessels, and parenchyma’s, the complex autocorrelation models based on different types of the dynamic scattering regimes $d$, amount of static scatterers $\rho$ and $\tau _c$ were generated. These were used for the generation of speckles using SDE and their auto-correlation was compared with the desired auto-correlation.
Fig. 4.
Fig. 4. Experimental validation of SDE-based simulation of laser speckles. (a) Illustration of the LSCI experimental system used for imaging mice brain. (b) The experimental and the simulated SDE-based speckle intensity were fitted for $e^{-t/\tau _c}$ model to obtain the decorrelation time $\tau _c$. (c) The simulated speckle intensity is integrated upto 1 ms exposure time of the camera. (d) The Experimental and simulated speckle contrast from the simulated speckle intensity at exposure time T=1.1 ms.
Fig. 5.
Fig. 5. Simulation of laser speckles for deep tissue blood flow. Schematic diagram of the DCS system simulated by the SDE algorithm, wherein a laser source (LD) is focused on the sample, and the detector captures the scattered intensity. The sample had the following optical and dynamic properties: $\mu _a=0.1 cm^{-1},\mu _s^{'}=0.12 cm^{-1} D_B=10^{-8} cm^2/s$. (b) The desired intensity has a diffusing profile from the centre (of the source) towards the boundary. (c) The simulated speckle intensity using the proposed method are shown at an exposure time of 1 ms. Both (b) and (c) show the log of the speckle intensity for better visualization. (d) The profile of the simulated speckles (in terms of $log(\phi.r)$) is fitted to a linear plot with the slope comparable to the desired slope indicating the match of optical properties. (e) The simulated speckle intensity integrated up to 10 $\mu s$ exposure time. (Refer to the supplementary video Visualization 3 for the intensity speckle fluctuations as a function of time and space). (f) The normalized intensity auto-correlation $g_2$ at two representative SD pairs, r =1 cm and 3 cm, along with the desired auto-correlation obtained using the semi-infinite solution of CDE. (g) The histogram of the simulated intensity speckle fits to a gamma distribution with the desired $\mu$. (h) and (i) shows the auto-correlation $g_2$ at delay times $\tau =10 \mu s$ and $100 \mu s$ as a function of spatial coordinates. (j) Speckle contrast is plotted for multiple exposure times at two representative SD pairs (r=1 cm and 3 cm).
Fig. 6.
Fig. 6. In-vivo animal and human experimental results for validating the proposed algorithm. (a) The autocorrelation $g_2$ of the experimental and simulated diffuse speckle intensity measured at an SD separation of 0.4cm in the brain of a mouse in alive and sacrificed condition. (b) The DCS measurement was performed in an adult human brain at an SD of 2 cm, and the normalized field auto-correlation was extracted based on Siegert’s relation. Intensity autocorrelation obtained from the simulated speckle is comparable to the experimental results.
Fig. 7.
Fig. 7. Effect of noise added to speckle intensity for realistic modelling of speckle intensity. (a) Intensity profile $log(\phi.r)$) is plotted as a function of $r$ at an exposure time of $1\mu s$. When no noise is added, the SDE generates an intensity profile that follows a linear fit. When the camera noise was added, the SNR decreases as the SD separation increases and the signal deviates from the linear fit. (b) At a given representative SD separation, speckle intensity generated by SDE follows the desired autocorrelation when no noise is added. The detector noise is applied (based on APD), resulting in the deviation of the autocorrelation from the desired value for larger values of delay times. (c) The simulated speckle contrast at various exposure times $T$ is plotted. At larger time $T$, due to saturation of the camera, the speckle contrast deviates from the desired value. (d) The mean and the variance of the multi-exposure speckle contrast calculated by varying the number of speckles $N$. As the number of speckles used to calculate speckle contrast increases, the standard deviation decreases and mean speckle contrast approaches the analytical solution.

Tables (1)

Tables Icon

Algorithm 1. SDE - Milstein scheme

Equations (6)

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

d x ( t ) = a ( x ( t ) , t ) d t + b ( x ( t ) , t ) d W ( t ) ; x ( 0 ) = x 0 .
( a p ) x = 1 2 2 ( b 2 p ) x 2 .
b 2 ( x ( t ) ) = 2 p ( x ) 0 x a ( x ( s ) ) p ( s ) d s = 2 p ( x ) 0 x α ( x ( s ) μ ) p ( s ) d s = 2 α μ x .
d I ( t ) = α ( I ( t ) μ ) d t + 2 α μ I ( t ) d W ( t ) ,
I n + 1 = I n + a ( I n ) Δ t + b ( I n ) Δ W n + 0.5 b ( I n ) b ( I n ) [ ( Δ W n ) 2 Δ t ] ,
g 2 ( τ ) = 1 + β ρ 2 ( d | g 1 ( n = X ) ( τ ) | + ( 1 d ) | g 1 ( n = 1 ) ( τ ) | ) 2 + β ( 1 ρ ) ρ ( d | g 1 ( n = X ) ( τ ) | + ( 1 d ) | g 1 ( n = 1 ) ( τ ) | ) + C
Select as filters


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