Recovery of the diffuse correlation spectroscopy data-type from speckle contrast measurements: towards low-cost, deep-tissue blood flow measurements.

A multi-step Volterra integral equation-based algorithm was developed to measure the electric field auto-correlation function from multi-exposure speckle contrast data. This enabled us to derive an estimate of the full diffuse correlation spectroscopy data-type from a low-cost, camera-based system. This method is equally applicable for both single and multiple scattering field auto-correlation models. The feasibility of the system and method was verified using simulation studies, tissue mimicking phantoms and subsequently in in vivo experiments.


Introduction
Quantitative in vivo imaging of blood flow using laser speckles have been investigated by several researchers in the past. Some of the techniques for quantifying blood flow includes, laser speckle contrast imaging (LSCI) [1][2][3], laser doppler flowmetry (LDF) [4,5], diffuse correlation spectroscopy (DCS) [6][7][8] and its three dimensional tomographic extension called diffuse correlation tomography (DCT) [9][10][11]. While LSCI and LDF have been used to quantify blood flow in relatively superficial layers of tissue, the deep tissue blood flow is quantified using DCS and DCT.
The seminal work to employ laser speckles for imaging blood flow was reported in Ref [12], where the time integrated speckles recorded using camera was used to photograph blood flow in human retina. Although limited to the imaging depth of approximately 1 mm, LSCI provided inexpensive and faster way of real time blood flow imaging with relatively simpler optical instrumentation. Based on LSCI method, multi-exposure recording of the speckle data was utilized for a quantitative recovery of absolute flow which is termed as multi-exposure speckle contrast imaging (MESI) [13]. In MESI, the measured multi-exposure speckle contrast images were fitted against single scattering flow models to quantify the flow.
Meanwhile the possibility of extending the dynamic light scattering to account for multiple scattering was explored in [14,15], which resulted in the development of diffusing wave spectroscopy (DWS). In DWS, under the ergodicity assumption, the temporal auto-correlation of the intensity speckles were utilized to measure the dynamics of the media [16,17]. Subsequently, the development of the correlation diffusion model led to diffuse correlation spectroscopy (DCS) which is extensively used in measuring deep tissue (> 1 cm) blood flow [8,[18][19][20].
A single channel DCS system comprises of a sensitive and expensive photo counting device like photon multiplier tube (PMT) or avalanche photo diode (APD) along with a hardware or software auto-correlator [6,8]. In this regard, usually multi channel acquisition systems are desired, as employing many detectors at a given spatial location improves the signal-to-noise ratio (SNR) of the system [21]. Also inorder to spatially resolve the blood flow, leading to a diffuse correlation tomography (DCT) system, many source detector (SD) pairs have been employed. For instance, in Ref. [9], 24 SD pairs were used to image over an area of 72 mm 2 , while Ref. [11] used 48 SD pairs to scan the area of 100 mm 2 in rat brain. One of the primary requirements of DCS/DCT system is many source detectors pairs which makes it a relatively expensive imaging modality. A solution to the above mentioned problem was proposed by Speckle Contrast Optical Spectroscopy (SCOS) by employing a camera, wherein each pixel of the camera acts a detector to simultaneously measure time integrated speckles for deep tissue flow imaging in Ref. [22,23], for both phantoms as well as invivo studies. SCOS employed a quantitative recovery of the flow by using the correlation diffusion model [6,8] and noise corrected speckle contrast measurements.
SCOS was further extended to achieve three dimensional tomographic imaging of blood flow as presented in Ref. [24] termed as speckle contrast optical tomography (SCOT). A high density version of SCOT was successfully applied to generate the tomographic images of the reduction in cerebral blood flow during local ischemic stroke in mice brain [25]. In Ref. [26], an electron multiplying charge coupled device (EMCCD) camera was employed to generate three dimensional flow images, in phantoms. Several successive publications have discussed the utility of employing array detectors like CCD or CMOS (complementary metal-oxide semiconductor) camera to probe tissue blood flow by measuring the speckle contrast data [27][28][29][30][31][32]. A compact version of SCOS using single photon avalanche diode (SPAD) array sensor was introduced in [33], wherein the flow was quantified using a multi-exposure speckle contrast measurement.
In a theoretical perspective, all the aforesaid methods were based on a relation connecting speckle contrast to normalized field auto-correlation [34] which in turn is related to flow. The speckle contrast is a time average over exposure of the normalized field auto-correlation. Hence, a multi-exposure speckle contrast data was used to quantify flow using least square minimization, without attempting to recover the field auto-correlation function [13,23,33]. A direct intensity auto-correlation measurement of blood flow is not possible due to the relatively larger exposure time of the camera which results in the time averaged measurements. Utilizing CCD / CMOS camera to directly measure intensity auto-correlation function has been attempted for DWS studies in Ref [35]. However the dynamics of the media utilized in this study, is slower when compared to dynamics of blood flow. Recently interferometric DCS has been attempted in Ref. [36], to measure blood flow using a CMOS camera, which is based on the principle of multi-mode fiber based interferometry.
In this paper, using a new algorithm, we present an inexpensive system using low frame rate CCD or CMOS cameras which can be employed for the recovery of the DCS data-type for invivo blood flow measurements. To the best of our knowledge, an inexpensive, low-frame rate camera has not been previously used to recover the full auto-correlation function, i.e. the DCS data-type. The algorithm is based on multi-step Volterra integral method (MVIM) [37] which uses multi-exposure speckle contrast measurement to recover the normalized field auto-correlation. Though the speckle contrast, at a given exposure, has averaged out the correlation decay, we retrieve the field auto-correlation function from multi-exposure speckle contrast data using MVIM. We demonstrate the system and method in imaging flow in tissue mimicking phantoms and also during human hand-cuff occlusion experiment. By recovering field auto-correlation function from multi-exposure speckle contrast data, we also prove that, speckle contrast infact contains information on the field auto-correlation. This fact is demonstrated by multi-exposure speckle contrast methods, as reported in previous publications [13,23,33], wherein the flow is quantified directly by least square fitting against speckle contrast.

Multi-step Volterra integral method (MVIM): theory and algorithm
The relation connecting speckle contrast (κ(r, T)) and normalized electric field auto-correlation function (g 1 (r, τ)) is given by Here, κ(r, T) is the speckle contrast given at a source detector separation r and at exposure time T, τ is the correlation delay time and β is the experimental constant accounting for the collection optics. In order to recover g 1 (r, τ) at each source detector separation r and for all relevant τ's, we pose the Eq. (1) in terms of volterra integral equation of first kind [38] as follows, where the term Ψ(T, τ) is called the kernel defined as Ψ(T, τ) ≡ (T − τ). With these definitions, the problem to recover field auto-correlation function using MVIM is stated as follows: Given κ 2 (r, T) for all T ∈ [T min , T max ] and for every source detector separation r, find g 1 (r, τ) for all τ ∈ [τ min , τ max ]. The selection of range of T's([T min , T max ]) and τ's([τ min , τ max ]) will be discussed later.
The standard method to solve volterra integral equations requires the kernel to be non-zero at T = τ [39]. Since the kernel in Eq. (2) violates this condition, we adopt the method given in Ref. [37] to solve the integral equation. In simpler terms, the method suggests to introduce a small shift in τ by δτ, so that T τ as they differ by δτ. The Eq. (2) is valid for every source detector (SD) separation r and hence we can drop the argument r from the function henceforth.
The MVIM algorithm to recover field auto-correlation was implemented numerically by discretizing the integral in Eq. (2) using trapezoidal integral rule as given below, Here, h is the step size, N is the number of discretization of the interval [0, T] and n is the index of the discretization. In matrix form, the above equation can be represented as, For simulation studies, speckle noise was added by using the noise model given in [32,40] . The standard deviation of the above-said noise model is given by σ κ = κ( where P is the number of pixels used for calculating speckle contrast. The matrix A is a lower triangular matrix and hence the best way to solve the above matrix equation is to apply the method of successive substitution. But due to the practical difficulties of measuring a dense multi-exposure speckle contrast data and the presence of associated noise, we adopted Tikhonov regularized least square minimization to solve the above matrix equation. The cost function ||Ax − b|| 2 2 + λ||x − x 0 || 2 2 was minimized for x, where λ is the regularization parameter determined by L-curve method [41] and x 0 is the prior information on g 1 (τ). Three parameters that play a key role in the above minimization problem are correlation delay time (τ), exposure time (T) and the prior information (x 0 ). The optimal selection of the above parameters for the better recovery of g 1 (τ), is explained in the following sections. Let the range of τ s at which we seek the normalized field auto-correlation function g 1 (τ), be [τ min : ∆τ : τ max ]. By looking at the limits of integration of Eq. (1), it is clear that τ should be restricted to T max , which is the maximum exposure time. Hence, the maximum allowed correlation delay time τ max = T max . The value of τ min has to be selected depending on the characteristic correlation delay time (τ c ) of the sample under measurement. To capture correlation decay, τ min should be selected such that τ min <<τ c .
The correlation delay time step ∆τ was determined by comparing the accuracy of the numerical solution of Eq. (1) against the analytical solution. The numerical speckle contrast κ N (T) was obtained by applying a trapezoidal integration to Eq. (1). The expression for analytical speckle contrast to Eq. (1) for single scattering case is given by Here the field auto-correlation model used was g 1 (τ) = e −τ/τ c [34]. We select ∆τ such that the percentage error between κ N and κ A was less than a pre-determined value.

Exposure time T
The range of the exposure time needed to generate multi-exposure speckle contrast data are often limited by the dynamic range of the camera. We adopted the speckle contrast sensitivity analysis, wherein the sensitivity is defined as S a = −T dκ d(T/τ c ) as given in [40,42] , to find an optimal exposure time. In Ref. [42], it was proved that the speckle contrast has maximum sensitivity when T equals τ c , the characteristic decay time of the sample. Although the above scenario was shown to hold for single scattering case as in LSCI, we have found that it is equally applicable for diffusing regime as well (see section 4.1). Hence T min was to be chosen as T min ≤ τ c , while T max was chosen to cover the whole correlation decay. We found that T max of approximately 100 time τ c (where g 1 (τ) → 0) serves the purpose.

Prior information
While minimizing the cost function to recover the field auto-correlation, the selection of prior, i.e., x 0 , plays a crucial role. When τ<<τ c , the value of normalized field auto-correlation function will tend to unity, i.e. g 1 (τ<<τ c ) ≈ 1 and hence we chose x 0 = 1 for appropriately small τ values. Note that the prior x 0 is a function of τ and hence by x 0 = 1 implies it has a value of unity for all τ's. This assumption, can lead to certain errors between the actual and reconstructed g 2 1 (τ) and hence we also adopt an iterative scheme to update the prior at each iteration.

Iterative algorithm for MVIM
As explained above, we adopt an iterative algorithm to update the prior x 0 at each iteration by fitting the recovered g 1 (τ) against the appropriate model. The flow chart of the iterative scheme is shown in Fig. 1. The algorithm starts by recovering g 1 (τ) using MVIM with a prior of x 0 = 1 for all relevant τ s. The recovered g 1 (τ) was fitted for D B or τ c for the cases of multiple scattering or single scattering respectively. The Correlation Diffusion Equation (CDE) was used as the model for multiple scattering [6] while single scattering utilizes simple model like g 1 (τ) = e −(τ/τ c ) [34]. The g f 1 corresponding to fitted D B or τ c was used as the prior to next iteration. This process was continued till the percentage error between two consecutive recovered g 1 (τ) lies below a tolerance level. Note that the prior information from second iteration will not be unity for all τ values, instead it will move closer to actual g 1 (τ) to be recovered. In case the model is unknown (as in the case of two layer or multi layer models), we may choose to fit for an empirical model:

Proposed system
The experimental setup to validate the proposed method to measure field auto-correlation function using an inexpensive low frame rate camera is given in Fig. 2. A current and temperature controlled laser diode (785 nm, 90 mW, Thorlabs), with beam shaping optics (aspheric lens, anamorphic prism, aperture and focusing lens) was used to form a pointed source of diameter less than 1 mm diameter. The beam was focused to the sample (tissue or phantoms) and the scattered intensity was measured using a CCD camera (Basler acA-640-120-um) in the reflection geometry. An objective lens of focal length of 50 mm and f-number f /# of 8 was used to match the speckle size to pixel size of the camera. Additionally, a LSCI system was developed for single scattering case wherein uniform illumination was produced by using diffusers.
The scattered intensity from the phantom or tissue at different exposure times (T) were recorded by the camera. The speckle contrast was calculated temporally over 500 frames and was corrected for dark and shot noise as given in [23,24]. To calculate speckle contrast κ(r, T) at a given Source Detector SD separation (r) with high SNR (Signal to Noise Ratio), we defined detectors in form of annuals (or rings) with radius r from the source with inner and outer diameters being r − 0.01cm and r + 0.01cm respectively [23]. The speckle contrast that falls within these detectors for a given r was averaged and was used as speckle contrast for a given SD separation r and this was repeated for every exposure time T.

Tissue mimicking phantoms
To create tissue mimicking phantoms of varying dynamic properties, two phantoms based on Intralipid (Fresenius Kabi India, Intralipid 20%) were made. The phantoms, one with distilled water as base (hereafter denoted as "Intralipid phantom") and another with glycerol as base (hereafter denoted as "Glycerol phantom") were prepared as per the recipe given in [23,43,44]. The absorption coefficient (µ a ) and the reduced scattering coefficient (µ s ) of both the phantoms were estimated by fitting against the analytical solution of diffusion equation [45]. The µ a and µ s of intralipid phantom (D a B ) was estimated to be 0.021 cm −1 and 8.1 cm −1 respectively and while that of glycerol phantom (D b B ) was 0.021 cm −1 and 10.5 cm −1 respectively. The value of β was determined by calculating κ at small exposure time T, such that T<<τ c . In addition, flow phantoms were made, where the flow was controlled using a syringe pump (NewEra pump systems Inc.). The liquid phantoms were made to flow on a semi circular shaped tube, embedded on silicon elastomer cavity (Slygard 184 mixed with TiO 2 ), with a diameter such that, the light does not interact with the boundaries of the supporting cavity and hence the semi-infinite geometry for light propagation can be adopted.

In vivo experiment
To show the feasibility of working of the proposed method to measure field auto-correlation function for in vivo tissues, we have performed a human forearm cuffing experiment, wherein the blood flow in human arm was occluded by using a standard sphygmomanometer. This work was undertaken with approval of the Institute Ethical Committee of the Indian Institute of Technology -Bombay, Approval Number: IITB-IEC/2018/017. The optical properties such as µ a and µ s were calculated from the scattered intensity data by fitting it against the solution of diffusion equation [46]. The multi-exposure speckle contrast data was acquired at pre-occlusion, occlusion and post-occlusion phases (measured after 120 s from releasing the cuff of sphygmomanometer). This multi-exposure speckle contrast data was used for recovering normalized field auto-correlation function (g 1 (τ)) using MVIM. The blood flow index (BFI), which is αD B [8], was quantified from the recovered field auto-correlation, where the quantification can be carried out using CDE model (Green's function solution used in this work) [8] or using other approximations such as the modified Beer Lambert law for DCS [47,48] or high order linear algorithms based on Monte Carlo using Taylor polynomial expansion [49,50].

Results
We have validated the MVIM algorithm, to measure field auto-correlation function using the camera based system, by both numerical simulations as well as experimental measurements. We have shown the numerical results for both single scattering and multiple scattering cases. Subsequently, we have also experimentally validated our method by tissue mimicking phantoms and human in vivo experiments.

Validation of MVIM using numerical simulations
In this section, we present the numerical simulations to validate MVIM for both single and multiple scattering regimes. The single scattering model for field auto-correlation utilizes g 1 (τ) = e −(τ/τ c ) as in LSCI [13,34], whereas, the CDE was used to generate field auto-correlation for multiple scattering as adopted in DCS [6]. Trapezoidal integration rule was adopted to evaluate the integral in Eq. (1) to compute the speckle contrast for the assumed auto-correlation models. As mentioned in Section 2, speckle noise model given in [32,40] was used to add noise to the simulated speckle contrast data.
The simulated multi-exposure speckle contrast data generated as mentioned above was fed to the MVIM algorithm presented in Section 2. We have chosen an optimal exposure time by computing the sensitivity of speckle contrast to exposure time. The sensitivity shows a maximum, when the exposure time equals the characteristic decay time, i.e., T = τ c [42]. Hence the multi-exposure speckle contrast data will be sensitive to the characteristic decay time of the sample, if the exposure time is selected in the vicinity of τ c . Here we have selected T such that, τ c ≤ T ≤ T max and utilize the speckle contrast data to recover g 1 (τ) using MVIM.
The range of correlation delay time τ, at which we seek the field auto-correlation function g 1 (τ), was selected to be [10 −10 : T max ]. The correlation delay time step ∆τ was determined by comparing the accuracy of the numerical solution (κ N ) of Eq. (1) against the analytical solution (κ N ) as mentioned in section 2. A plot of percentage error (E) between the κ N and κ A for a single scattering case against the number of τ's (in the range of 10 −10 to 1 ms with τ c = 1 ms) is shown in Fig. 3. We have used N = 250, where the percentage error E is almost negligible, as shown in Fig. 3.   Fig. 3. Plot between Number of τ s, N, and percentage error (E = κ A −κ N κ A ×100) indicating that at around 250 intervals the percentage error is almost negligible (less than 0.002%) The dependence of recovered g 1 (τ) on the prior information x 0 was explained in the section 2, where the rationale for adopting an iterative scheme was mentioned. We present simulation results of the MVIM algorithm under following special cases: (a) With and without prior i.e., x 0 = 1 and x 0 = 0 respectively, (b) with an optimal exposure range i.e., τ c ≤ T ≤ T max and (c) iterative scheme as explained in section 2 (see Fig. 1). In addition, we also present the recovery of g 1 (τ) for a two layer model (i.e. particle diffusion coefficient D B assumes two different values in two different layers) wherein we computed the corresponding speckle contrast data using a Finite Element Method (FEM) based forward solver for CDE [51]. Here g 1 (τ) exhibits two different decay characteristics corresponding to each layers which was also recovered successfully using MVIM.

Single scattering
In the single scattering case, we have used the normalized field correlation model of the form g 1 (τ) = e (−τ/τ c ) [13,34]. We have used two different τ c values, which are τ a c = 1 ms and τ b c = 100 ms to generate g 1 (τ) and κ(T) from the analytical models [34]. Note that we have added speckle noise only to the speckle data generated for τ b c = 100 ms to differentiate the cases with and without noise. For the samples, with characteristic time τ a c and τ b c , the exposure ranges, i.e. T min to T max , used were 1 µs to 1 s and 1 µs to 100 s respectively. The number of exposures used was 50, spaced equally in the above exposure range in logarithmic scale for both the samples. Figure 4(a) shows the recovered g 1 (τ) from κ(T) for two different τ c values, where the prior information is chosen to be x 0 = 0. Here the exposure time was not chosen optimally with respect to τ c . Instead, a wide range of T, such that T min ≤ T<τ c ≤ T max was used to generate the simulated multi-exposure speckle contrast data. The recovered g 1 (τ) shows a transient from 0 to unity for lower τ values (τ<T min , which was 10 −6 s) due to fact that x 0 is chosen to be zero. Figure 4(b) shows similar plots for the case, where x 0 = 1. The initial transient of g 1 (τ) as seen in Fig. 4(a) is absent here because of the fact that prior is now at x 0 = 1, which causes the minimization problem to search for an optimal solution in the vicinity of g 1 (τ) = 1. Measuring speckle contrast (κ(T)) for a wide range of exposure time is practically not feasible and hence we adopt the sensitivity analysis as described before to find an optimal exposure time. The sensitivity curves which plots S a against T/τ c is shown in Fig. 5(a). Here the range D shows the optimal set of exposure times T, where T ≥ τ c . The speckle contrast data is plotted against the exposure time in Fig. 5(b), wherein the range of optimal exposure time (i.e., range D in Fig. 5(a)) used to recover g 1 (τ) has been highlighted. By employing MVIM, we recover g 1 (τ) from the above mentioned optimal exposure time T, such that T ≥ τ c and the recovered g 1 (τ) is shown in Fig. 6(a). It shows several ripples specifically for larger τ values, which is due to the fact that number of τ's at which the g 1 (τ) is sought, is larger than the number of exposure values used. Note that, in this case, the prior information x 0 = 1 was used in recovery of g 1 (τ). The noise added to multi-exposure data and the prior information (x 0 = 1) also contributes to the above mentioned ripples.
In order to minimize the ripples, we decided to adopt an iterative algorithm for MVIM as discussed in Section 2.2, where the prior was updated by fitting g 1 (τ) against single scattering model (see Fig. 1). The results of the iterative scheme are shown in Fig. 6(b), where amount of ripples is reduced. The residual between the original and the recovered g 1 (τ) is shown in Fig. 7(a) and 7(b) for the above-said two τ c values by using iterative scheme. The recovered g 1 (τ) is fitted for τ c using the single scattering model for field auto-correlation and the residual norm  between original and recovered g 1 (τ), is tabulated in Table 1. The results in the Table 1 suggest that employing MVIM using appropriate prior and optimal exposure along with an iterative scheme, ensures a better recovery of field auto-correlation function g 1 (τ) from multi-exposure speckle contrast data. Fig. 7. Residual between original and recovered g 1 (τ) using MVIM by iterative scheme (plot in Fig. 6(b)), is shown for (a) τ a c , i.e., noiseless data (b) τ b c , i.e., noisy data case

Multiple scattering
The analytical solution of CDE was used for simulating the normalized field auto-correlation function for multiple scattering case [6,46]. The µ a and µ s used were 0.1 cm −1 and 8 cm −1 respectively and SD separation of 1 cm was used. Two different particle diffusion coefficients (D B ), i.e. D a B = 2 X 10 −8 cm 2 /s and D b B = 2 X 10 −10 cm 2 /s were used and the speckle noise was Table 1. Table comparing original and recovered τ c values and residual norm (l 2 ), where the first two rows shows the recovered τ c without and with prior and the third row shows the recovered τ c , with speckle contrast data in optimal exposure range. The fourth row shows recovered τ c and residual norm for the iterative scheme.

Without noise
With noise added to the latter (speckle contrast data corresponding to D b B ) using the noise model mentioned in Section 2. The exposure range i.e., T min to T max used were 0.01 µs to 0.1 s and 1 µs to 100 s for D a B and D b B samples respectively with 50 exposures equally spaced in logarithmic scale. Figure 8 shows the results of MVIM method, where g 1 (τ) has been recovered from multiexposure speckle contrast data. Figure 8(a) shows the recovered normalized field auto-correlation using MVIM, when no prior information, i.e., x 0 = 0, was used for two different D B values and it can be seen that there is transient curve from 0 to 1 at τ ≤ T min . Different exposure ranges were used for different samples to demonstrate the initial transient curve during recovery, when no prior information, i.e., x 0 = 0 was used. Here, from Fig. 8(a), we see that g 1 (τ) can only be recovered for τ s, such that τ ≥ T min , when no prior information i.e. x 0 = 0 is used. Similarly, Fig. 8(b) shows the plot for the case where the prior information x 0 = 1 was used. As seen in Fig. 8(b), the initial transient curve is absent due to use of prior information x 0 = 1. Here, a wide range of exposure time T, i.e., T min ≤ T ≤ T max were used.
We adopt the sensitivity analysis as described in single scattering case, where S a is plotted against T/τ c as shown in Fig. 8(c). Although the relation between D B and τ c has been derived in [48], we choose to follow a simpler method of choosing τ c = τ such that g 1 (τ) = 1/e ≈ 0.37. The τ c values were 50 µs and 5 ms for the samples D a B and D b B respectively. The optimal exposure time T is chosen such that τ c ≤ T ≤ T max and g 1 (τ) is recovered as shown in Fig. 8(d). To utilize the prior information in an adaptive manner for better reconstruction, the iterative scheme as described in Section 2 ( Fig. 1) was used. The results of iterative scheme based recovery of g 1 (τ) is shown in Fig. 8(e). The residual between the original and recovered g 1 (τ) for the noisy speckle contrast data (iterative scheme) is shown in Fig. 8(f). The g 1 (τ) recovered was fitted for D B by using CDE model and is tabulated in Table 2, along with the residual norm between the original and recovered g 1 (τ). where the first two rows shows the recovered D B without and with prior and the third row shows the recovered D B , with speckle contrast data in optimal exposure range. The fourth row shows recovered D B for the iterative scheme.   [42], where τ c was chosen as time τ, when g 1 (τ) = 1/e ≈ 0.37. From this curve it can be seen that maximum sensitivity is achieved, when T ≈ τ c ; Recovered g 1 (τ) from κ(T) by (d) using optimal exposure range i.e., T ≥ τ c (e) iterative scheme; (f) Residual between original and recovered g 1 (τ) using iterative scheme. Here the superscript M in the legend indicates that it is MVIM recovered and O indicates that it is original g 1 (τ)

Two layer model
In order to validate our method in two layer tissue models, we have simulated the auto-correlation function, g 1 (τ) using FEM based forward solver [51] for CDE. The sample has two different D B values for two layers (first layer up to 1 cm thickness with D B = 2 × 10 −12 cm 2 /s and the other layer of D B = 2 × 10 −8 cm 2 /s) and the SD separation was 2 cm with µ a = 0.1 cm −1 and µ s = 8 cm −1 . The multi-exposure speckle contrast (κ(T)) was computed numerically using Eq. (1) and statistical noise was added [42]. The exposure range used was 10 −6 s to 10 2 s with 50 exposures spaced in logarithmic scale with a decorrelation time (τ c ) of 0.4 ms. Figure 9 shows the results of MVIM based recovery of g 1 (τ) for a two layer model. Figure 9(a) shows the original and recovered g 1 (τ) using MVIM, when prior information, x 0 = 0 was used. Figure 9(b) and (c) shows the similar plots, wherein g 1 (τ) was recovered from the noise added speckle contrast, by using prior information x 0 = 1 and iterative scheme respectively. Figure 9(d) shows the sensitivity curve, as discussed in single scattering case, with two maxima's indicating that there are two different flows (or D B values).

Single scattering
The experimental system shown in Fig. 2 (as explained in section 3) was modified to generate LSCI data (κ(T)). The speckle contrast was calculated from 500 images with 20 exposures in equally spaced in logarithmic scale with different exposure ranges for different flows due to limited dynamic range of the camera. For flow τ a c , the exposure range was 0.6 ms to 8 ms. Similarly for τ b c and τ c c , the exposure ranges used were, 1.5 ms to 16 ms and 5 ms to 100 ms respectively. From the multi-exposure speckle contrast data (κ(T)), by using MVIM, g 1 (τ) was recovered and was fitted for τ c using the analytical model of g 1 = exp(−τ/τ c ). Figure 10(a) shows the recovered g 1 (τ) using MVIM with prior information, x 0 = 1, and compares it with g 1 (τ), computed for τ c fitted by MESI [13]. Figure 10(b) shows similar plot, wherein iterative scheme was employed. As the flow increases, the g 1 (τ) shifts towards left indicating the increase in flow. Table 3 compares the τ c values of MESI and MVIM based methods and it can be seen that they are in reasonable agreement with one another.

Multiple scattering
A laser source was focused on to the sample to form a point source as shown in Fig. 2 and the multi-exposure speckle contrast data (κ(T)) was measured. As mentioned in Section 3, the speckle contrast was measured for a SD separation r of 0.5 and 1 cm for both intralipid and glycerol phantoms at different exposures. For the calculation of speckle contrast, we have used 500 frames with 20 exposure times ranging from 0.05 ms to 2 ms for intralipid phantom and 1 ms to 100 ms for glycerol phantom. The multi-exposure speckle contrast data as a function of exposure time at two different SD separations are shown in Fig. 11. Using MVIM, the normalized field auto-correlation function was recovered for both intralipid and glycerol phantoms for the above mentioned SD separations and the results are shown in Fig. 12 and Fig. 13. It is compared with g 1 (τ) computed for D B fitted by SCOS [23]. Figure 12 shows the recovered g 1 (τ) by MVIM using prior information x 0 = 1, where (a) shows the results of intralipid phantom and (b) shows the results of glycerol phantom. Similar plot is shown in Fig. 13 by using iterative scheme, wherein solution of CDE [6,46] was used as model for updating the prior information. Recovered g 1 (τ) using MVIM, using prior (x 0 = 1) and compared against g 1 (τ) fitted using SCOS with SD -0.5 and 1 cm for (a) intralipid phantom (b) glycerol phantom.
The recovered g 1 (τ) was fitted for D B against CDE model [6,46]. For Intralipid phantom, with prior information (x 0 = 1), the fitted D B values are 1.46 X 10 −8 cm 2 /s and 1.04 X 10 −8 cm 2 /s Subsequently, using flow phantoms controlled by syringe pump, different flows were made and speckle contrast at different exposures were measured at an SD separation of 0.75 cm. By using MVIM, normalized field auto-correlation (g 1 (τ)) function was recovered and is plotted in Fig. 14(a) for different flows. Here iterative MVIM algorithm was employed to recover g 1 (τ). The g 1 (τ) was fitted for D B values against solution of CDE and plotted against the flow velocity is shown in Fig. 14(b) where it can be seen that the D B value increases linearly with velocity of the syringe pump. The D B fitted for MVIM is compared against the D B fitted for SCOS and is tabulated in Table 4.

in vivo experiment
We validate the proposed system and method in vivo, by performing a human hand blood cuff experiment, for 5 healthy volunteers. The experimental protocol involves 90 s measurement of multi-exposure speckle contrast for each pre-cuff, cuff and post cuff phases. The multi-exposure Table 4

. Comparison between MVIM fitted D B and SCOS fitted D B for different flow rates. The superscripts 'M' and 'S' refers to the MVIM and SCOS based fit for D B values respectively, which
corresponds to the curves shown in Fig. 14(a).
data was acquired at 15 exposures in the range of 0.5 ms to 20 ms spaced in logarithmic scale, with 200 frames in each exposure to calculate the speckle contrast. The results of MVIM recovered g 1 (τ) is plotted, along with results of g 1 (τ) obtained by D B fitted for SCOS at a SD separation of 1 cm, in Fig. 15(a). The relative blood flow rBF [33] is plotted for each phase in Fig. 15(b) . It can be seen that on an average for 5 volunteers, there is almost 70 − 80% decrease in blood flow between pre-cuff and cuff phases and the flow is recovered back during post cuff phase, which is comparable with results reported in Ref. [23,33]. In nut-shell, we have validated the proposed MVIM algorithm to measure field auto-correlation function using a camera based system, for tissue mimicking phantoms and for in vivo experiments.

Discussion and conclusion
An algorithm and formalism for the recovery of the DCS data-type from speckle contrast measurements using a system based on an inexpensive CCD/CMOS camera to measure blood flow is presented. A direct auto-correlation on the measured intensity speckles requires a camera with high frame rate (typically few MHz) along with a good sensitivity (Quantum Efficiency (QE) > 50% and dynamic range > 30000 e − ), high SNR (readout noise < 2e − ) and a wide range of exposure control (typically 250 ns to few seconds). Although high frame rate cameras are commercially available (with typically around 10 6 FPS), there is always a trade-off between frame rate and SNR. However, in Ref. [52], the intensity correlation was measured for laser speckle contrast imaging using camera of high frame rate (20000 Hz). The measured field auto correlation was used to understand the appropriate model for the blood flow. This is one of the potential applications where a direct recovery of the autocorrelation function is required instead of a least square fitted blood flow index (BFI). In the above method, the authors have employed a high frame rate of 20000 fps due to the fact that the de-correlation time (τ c ) is in the order of milliseconds. By recovering DCS data-type using MVIM method and an inexpensive low frame rate camera, we circumvent the need of a high frame rate camera. The MVIM method utilizes the fact that speckle contrast is a definite integral of normalized field auto-correlation function over correlation delay time integrated upto exposure time of the camera. We identified this integral equation as Volterra integral equation of first kind and hence we adopted a multi-step method to solve it. The proposed multi-step volterra integral method (MVIM) was successfully integrated with the camera based laser imaging system to recover the field auto-correlation from the multi-exposure speckle contrast data. The camera employed was a CCD camera with a low frame rate (<120 frames), but with reasonably good sensitivity (QE of 42%, dynamic range 17800e − ) and exposure control (4µs to 10 s).
This method is equally applicable for measuring superficial tissue blood flow (LSCI) as well as deep blood flow (DCS), as the relation connecting speckle contrast and field auto-correlation is common for both the cases. Although the single exposure speckle contrast integrates the field auto-correlation over correlation delay time (upto the exposure time of camera), the multiexposure speckle contrast retains the details of auto-correlation decay of the sample. This fact is utilized by MVIM for recovery of field auto-correlation function g 1 (τ). Methods such as MESI [13] for single scattering case and SCOS [22,23,33] for multiple scattering case utilize multi-exposure speckle contrast data. These methods have quantified blood flow either through τ c or D B by a direct least square fit on the multi-exposure speckle contrast data without attempting to recover the field auto-correlation. Infact, by using the proposed MVIM based system, we have proved that a direct recovery of g 1 (τ) for every τ values is possible instead of quantifying flow by a least square fit on multi-exposure speckle contrast data.
The volttera integral equation connecting κ(T) and g 1 (τ) was discretized numerically and posed as matrix equation. Although the volterra integral equation of first kind is not an ill-posed problem [53], due to the presence of noise in the measured speckle contrast and practical difficulties of measuring κ(T) for dense and wide range of exposure times, Tikhonov regularization based non-linear least square minimization was used for solving the above equation. The optimal measurement and minimization parameters such as optimal exposure time and sampling interval of correlation delay time has been presented along with an iterative scheme for better recovery of field auto-correlation function. The simulation and experimental results, that includes recovery of g 1 (τ) in flow phantoms and in vivo experiments, validates the presented MVIM based system to measure field auto-correlation function.
The results presented in this paper shows that our method is infact comparable to existing speckle contrast methods. There are two aspects of the presented work that we would like to mention, which are "inexpensive" nature of the measurement system compared to conventional DCS and preliminary study on the "information content" of the multi exposure speckle contrast data. The former is evident from the fact that by using the speckle contrast as the measurand, we do not need a fast and sensitive detector array, and, yet, still recover the full DCS data-type. For addressing the second aspect of the proposed method namely 'information content' of speckle contrast data, we have proved via MVIM that a multi exposure data has enough information to retrieve the DCS data-type for a relevant range of delay times (at each source detector SD separation).
Though we have shown the MVIM method for recovering the normalized field auto-correlation function g 1 (τ), we can also recover the normalized intensity auto-correlation, g 2 (τ), by posing the integral equation as (κ 2 (T) + 1) = 2 T ∫ T 0 (1 − τ T )g 2 (τ)dτ. Note that this equation is independent of β and any uncertainties in measuring β will not affect the recovery of g 2 (τ).
We have experimentally validated the working of the proposed method using SD separations of 0.5 cm to 1 cm which are shorter when compared to the usually used SD separations in DCS for measurements on the adults (i.e., 2cm to 3 cm). This is due to the current limitations of the camera. However, MVIM method is applicable for multi exposure speckle contrast data at every SD separation. We may use a better camera with a higher sensitivity (in terms of quantum efficiency and full well capacity) to achieve a higher SD separations, still staying on the in-expensive side. We note that these SD separations are relevant in terms of animal imaging studies [10,11,25]. We have demonstrated the working of two layer models via simulations, however, experimental validation was not performed. This is primarily due to the low dynamic range of the camera used for the experiment. Since, two different correlation decay demands the data acquisition for a larger range of exposure time, the dynamic range acts as a constraining factor. We also note that the proposed non-iterative MVIM method does not require any field auto-correlation model (such as CDE or e −τ/τ c ) for recovering field auto-correlation function g 1 (τ). We have adopted the above mentioned models in this paper for validating the results against MESI and SCOS.
In order to acquire the multi-exposure speckle contrast data, we need to frequently vary the exposure time. This can be achieved by either (a) keeping the camera exposure time constant but limiting the laser exposure using an external modulator [13] or (b) by changing the camera exposure time. Here, we have used the second approach which is not tedious since most modern cameras allow this to be controlled in an automated manner by software. One of the limitations of the system is, that the low frame rate of the camera constrains performing certain experiments, where in the flow parameters to be measured changes rapidly with time. For instance, we were not able to retrieve the reactive hyperaemia, during blood cuff experiment as we need to measure multi-exposure speckle contrast data within a short duration of time. The current implementation of calculating speckle contrast is slow as we calculate the speckle contrast over time. Alternatively, it can be calculated over space so that only one frame is needed per exposure time. With the current frame-rate, this would allow data acquisition at a much faster rate (approximately 2 seconds). One of the limitations of the MVIM method is that, even though we use iterative scheme, we could not completely eliminate the ripples present in the recovered field auto-correlation at the larger correlation delay time (τ values). This can be improved by averaging more frames (as the SNR increases by square root of number of frames) and by using a high sensitive camera. The above-mentioned limitations of both system and the method have to be improved in order to adapt the system for clinical studies.