Speckle contrast optical tomography: A new method for deep tissue three-dimensional tomography of blood flow

A novel tomographic method based on the laser speckle contrast, speckle contrast optical tomography (SCOT) is introduced that allows us to reconstruct three dimensional distribution of blood flow in deep tissues. This method is analogous to the diffuse optical tomography (DOT) but for deep tissue blood flow. We develop a reconstruction algorithm based on first Born approximation to generate three dimensional distribution of flow using the experimental data obtained from tissue simulating phantoms. © 2014 Optical Society of America OCIS codes: (170.3880) Medical and biological imaging; (110.3010) Image reconstruction techniques; (110.6150) Speckle imaging; (110.6955) Tomographic imaging. References and links 1. A. Devor, S. Sakadžić, V. Srinivasan, M. Yaseen, K. Nizar, P. Saisan, P. Tian, A. Dale, S. Vinogradov, M. Franceschini, and D. A. Boas, “Frontiers in optical imaging of cerebral blood flow and metabolism,” J. Cereb. Blood Flow Metab.32, 1259–1276 (2012). 2. M. J. Leahy, J. G. Enfield, N. T. Clancy, J. ODoherty, P. McNamara, and G. E. Nilsson, “Biophotonic methods in microcirculation imaging,” Med. Laser Appl. 22, 105–126 (2007). 3. T. Durduran, R. Choe, W. B. Baker, and A. G. Yodh, “Diffuse optics for tissue monitoring and tomography,” Rep. Prog. Phys. 73, 076701 (2010). 4. T. Durduran and A. G. Yodh, “Diffuse correlation spectroscopy for non-invasive, micro-vascular cerebral blood flow measurement,” NeuroImage 85, 51–63 (2014). 5. C. Riva, B. Ross, and G. B. Benedek, “Laser doppler measurements of blood flow in capillary tubes and retinal arteries,” Invest. Ophthalmol. Visual Sci. 11, 936–944 (1972). 6. M. Stern, “In vivo evaluation of microcirculation by coherent light scattering.” Nature 254, 56–58 (1975). 7. V. Rajan, B. Varghese, T. G. van Leeuwen, and W. Steenbergen, “Review of methodological developments in laser doppler flowmetry,” Lasers Med. Sci. 24, 269–283 (2009). 8. A. F. Fercher and J. D. Briers, “Flow visualization by means of single-exposure speckle photography,” Opt. Commun.37, 326–330 (1981). 9. A. K. Dunn, “Laser speckle contrast imaging of cerebral blood flow,” Ann. Biomed. Eng. 40, 367–377 (2012). 10. D. A. Boas, L. E. Campbell, and A. G. Yodh, “Scattering and imaging with diffusing temporal field correlations,” Phys. Rev. Lett. 75, 1855–1858 (1995). 11. C. Zhou, G. Yu, D. Furuya, J. Greenberg, A. J. Yodh, and T. Durduran, “Diffuse optical correlation tomography of cerebral blood flow during cortical spreading depression in rat brain,” Opt. Express 14, 1125–1144 (2006). #205756 $15.00 USD Received 31 Jan 2014; revised 10 Mar 2014; accepted 11 Mar 2014; published 28 Mar 2014 (C) 2014 OSA1 April 2014 | Vol. 5, No. 4 | DOI:10.1364/BOE.5.001275 | BIOMEDICAL OPTICS EXPRESS 1275 12. D. A. Boas, “Diffuse photon probes of structural and dynamical properties of turbid media: theory and biomedical applications,” Ph.D. thesis, University of Pennsylvania (1996). 13. H. M. Varma, A. K. Nandakumaran, and R. M. Vasu, “Study of turbid media with light: Recovery of mechanical and optical properties from boundary measurement of intensity autocorrelation of light,” J. Opt. Soc. Am. A 26, 1472–1483 (2009). 14. H. M. Varma, B. Banerjee, D. Roy, A. K. Nandakumaran, and R. M. Vasu, “Convergence analysis of the newton algorithm and a pseudo-time marching scheme for diffuse correlation tomography,” J. Opt. Soc. Am. A 27, 259– 267 (2010). 15. N. Hyvönen, A. K. Nandakumaran, H. M. Varma, and R. M. Vasu, “Generalized eigenvalue decomposition of the field autocorrelation in correlation diffusion of photons in turbid media,” Math. Meth. Appl. Sci. (2012). 16. J. P. Culver, T. Durduran, D. Furuya, C. Cheung, J. H. Greenberg, and A. G. Yodh, “Diffuse optical tomography of cerebral blood flow, oxygenation, and metabolism in rat during focal ischemia,” J. Cereb. Blood Flow Metab. 23, 911–924 (2003). 17. G. Dietsche, M. Ninck, C. Ortolf, J. Li, F. Jaillon, and T. Gisler, “Fiber-based multispeckle detection for timeresolved diffusing-wave spectroscopy: characterization and application to blood flow detection in deep tissue,” Appl. Opt.46, 8506–8514 (2007). 18. T. Binzoni, T. S. Leung, D. Boggett, and D. Delpy, “Non-invasive laser doppler perfusion measurements of large tissue volumes and human skeletal muscle blood rms velocity,” Phys. Med. Biol. 48, 2527 (2003). 19. S. R. Arridge and J. C. Schotland, “Optical tomography: forward and inverse problems,” Inverse Problems 25, 123010 (2009). 20. Y. Zhan, A. T. Eggebrecht, J. P. Culver, and H. Dehghani, “Image quality analysis of high-density diffuse optical tomography incorporating a subject-specific head model,” Front. Neuroenerg. 4, 103389 (2012). 21. V. Viasnoff, F. Lequeux, and D. J. Pine, “Multispeckle diffusing-wave spectroscopy: a tool to study slow relaxation and time-dependent dynamics,” Rev. Sci. Instrum. 73, 2336–2344 (2002). 22. A. P. Y. Wong and P. Wiltzius, “Dynamic light scattering with a ccd camera,” Rev. Sci. Instrum. 64, 2547–2549 (1993). 23. J. D. McKinney, M. A. Webster, K. J. Webb, and A. M. Weiner, “Characterization and imaging in optically scattering media by use of laser speckle and a variable-coherence source,” Opt. Lett. 25, 4–6 (2000). 24. B. J. Ackerson, R. L. Dougherty, N. M. Reguigui, and U. Nobbmann, “Correlation transferapplication of radiative transfer solution methods to photon correlation problems,” J. Thermophys. Heat Transf. 6, 577–588 (1992). 25. R. L. Dougherty, B. J. Ackerson, N. M. Reguigui, F. Dorri-Nowkoorani, and U. Nobbmann, “Correlation transfer: development and application,” J. Quant. Spectrosc. Radiat. Transfer 52, 713–727 (1994). 26. D. A. Boas and A. G. Yodh, “Spatially varying dynamical properties of turbid media probed with diffusing temporal light correlation,” J. Opt. Soc. Am. A 14, 192–215 (1997). 27. S. A. Carp, N. Roche-Labarbe, M.-A. Franceschini, V. J. Srinivasan, S. Sakadžić, and D. A. Boas, “Due to intravascular multiple sequential scattering, diffuse correlation spectroscopy of tissue primarily measures relative red blood cell motion within vessels,” Biomed. Opt. Express 2, 2047 (2011). 28. R. Bonner and R. Nossal, “Model for laser doppler measurements of blood flow in tissue,” Appl. Opt, 20, 2097– 2107 (1981). 29. R. C. Mesquita, T. Durduran, G. Yu, E. M. Buckley, M. N. Kim, C. Zhou, R. Choe, U. Sunar, and A. G. Yodh, “Direct measurement of tissue blood flow and metabolism with diffuse optics,” Philos. Trans. R. Soc., A 369 4390–4406 (2011). 30. S. R. Arridge, “Optical tomography in medical imaging,” Inverse Problems 15, R41 (1999). 31. C. Xu and Q. Zhu, “Light shadowing effect of large breast lesions imaged by optical tomography in reflection geometry,” J. Biomed. Opt. 15, 036003 (2010). 32. J. D. Briers, “Laser doppler, speckle and related techniques for blood perfusion mapping and imaging,” Physiol. Measurement 22, R35 (2001). 33. R. Bandyopadhyay, A. S. Gittings, S. S. Suh, P. K. Dixon, and D. J. Durian, “Speckle-visibility spectroscopy: A tool to study time-varying dynamics,” Rev. Sci. Instrum. 76, 093110 (2005). 34. M. v. van Rossum and T. M. Nieuwenhuizen, “Multiple scattering of classical waves: microscopy, mesoscopy, and diffusion,” Rev. Mod. Phys. 71, 313 (1999). 35. V. A. Markel and J. C. Schotland, “On the convergence of the born series in optical tomography with diffuse light,” Inverse Problems23, 1445 (2007). 36. B. C. White, “Developing high-density diffuse optical tomography for neuroimaging,” Ph.D. thesis, Washington University in St. Louis (2012). 37. J. P. Culver, A. M. Siegel, J. J. Stott, and D. A. Boas, “Volumetric diffuse optical tomography of brain activity,” Opt. Lett.28, 2061–2063 (2003). 38. S. Yuan, “Sensitivity, noise and quantitative model of laser speckle contrast imaging,” Ph.D. thesis, Tufts University (2008). 39. M. A. O’Leary, “Imaging with diffuse photon density waves,” Ph.D. thesis, University of Pennsylvania (1996). 40. H. He, Y. Tang, F. Zhou, J. Wang, Q. Luo, and P. Li, “Lateral laser speckle contrast analysis combined with line beam scanning illumination to improve the sampling depth of blood flow imaging,” Opt. Lett. 37, 3774–3776 #205756 $15.00 USD Received 31 Jan 2014; revised 10 Mar 2014; accepted 11 Mar 2014; published 28 Mar 2014 (C) 2014 OSA1 April 2014 | Vol. 5, No. 4 | DOI:10.1364/BOE.5.001275 | BIOMEDICAL OPTICS EXPRESS 1276 (2012). 41. J. F. Dunn, K. R. Forrester, L. Martin, J. Tulip, and R. C. Bray, “A transmissive laser speckle imaging technique for measuring deep tissue blood flow: an example application in finger joints,” Lasers Surg. Med. 43, 21–28 (2011). 42. A. Mazhar, D. J. Cuccia, T. B. Rice, S. A. Carp, A. J. Durkin, D. A. Boas, and B. J. T. B. Choi, “Laser speckle imaging in the spatial frequency domain,” Biomed. Opt. Express 2, 1553–1563 (2011). 43. R. Bi, J. Dong, and K. Lee, “Deep tissue flowmetry based on diffuse speckle contrast analysis,” Opt. Lett. 38, 1401–1403 (2013). 44. R. Bi, J. Dong, and K. Lee, “Multi-channel deep tissue flowmetry based on temporal diffuse speckle contrast analysis,” Opt. Express 21, 22854–22861 (2013). 45. M. Heckmeier, S. E. Skipetrov, G. Maret, and R. Maynard, “Imaging of dynamic heterogeneities in multiplescattering media,” J. Opt. Soc. Am. A 14, 185–191 (1997).


Introduction
Over the years, noninvasive optical imaging of microvascular blood flow has received much attention in biomedical research [1][2][3][4]. Laser Doppler flowmetry (LDF) [5][6][7] and laser speckle flowmetry (LSF) [8,9] are two promising techniques for two-dimensional, relatively superficial imaging of blood flow using dynamic laser speckles. Diffuse correlation spectroscopy (DCS) [3,10] based on similar physical principles as LSF and LDF but for deep tissues has been extended to three-dimensional tomography (diffuse correlation tomography (DCT) [11][12][13][14][15][16]). However DCT instrumentation is limited by fairly low signal-to-noise, low dynamic range in terms of detectable intensity and relatively expensive detectors that have kept the detector channel counts to less than ∼28 [17]. In this work, we present a new speckle contrast based tomographic approach-speckle contrast optical tomography (SCOT)-that provides efficient three-dimensional imaging of flow in turbid media that can be expanded to significantly large detector channel counts.
In LDF, the presence of moving scattering centers in a turbid medium, such as human tissues, causes a Doppler shift in frequency of the re-emitted light. The optical mixing of the Dopplershifted and non-shifted frequencies at the detector gives rise to a randomly fluctuating intensity distribution whose power spectral density is related to the flow of scatterers. For all the methods we discuss in this paper, the most significant contribution to the signal comes from the motion of red blood cells. For a typical separation of ∼250-1000 µm between the transmitting and receiving optical fibers the maximum penetration depth that can be achieved by LDF is ∼ 1 mm, thus limiting the technique to the superficial layers of tissue. There have been some proposals to utilize larger separation LDF probes based on multiple scattering models of quasielastic light scattering to extend the measurement depth of LDF to about ∼5 mm beneath the tissue surface, but they are not widely utilized [18].
In LSF, a camera (CCD or CMOS), records the scattered light emanating from a tissue that is uniformly or broadly illuminated with a laser source. The presence of moving scatterers induces time varying random fluctuations of the scattered intensity resulting in dynamic speckles. The camera integrates the dynamic speckle pattern during the exposure time and gives an integrated spatial intensity distribution. The spatial and temporal blurring of this integrated intensity pattern is studied using a statistical quantity called speckle contrast. Since LSF uses a uniformly illuminated laser source for probing, the penetration depth is limited to superficial layers of the tissue (typically less than 1 mm).
DCS, on the other hand, uses point sources and a photon diffusion model to probe deep tissues. DCT is normally carried out in an analogous fashion to diffuse optical tomography (DOT) [3,19]. It is a model based, computed tomography method where advanced inverse problem methods are utilized. Since, DOT has emerged as a promising technology in a variety of fields ranging from optical mammography to functional neuroimaging where it primarily measures parameters such as blood oxygen saturation and blood volume as well as the distribu-tion of contrast agents, we expect that a practical 3D tomography of blood flow in large tissue volumes would find rapid acceptance in the field.
The main disadvantage of DCS, and hence DCT, is the low signal-to-noise ratio (SNR) and the limited dynamic range which arises due to the fact that the intensity autocorrelation function has to be computed by measuring each speckle independently using single mode or few-mode fibers with small collection areas, i.e. collecting photons of order 10,000 per second. Since SNR increases as the square root of number of speckles, obtaining significant benefits in SNR as well as the dynamic range by doing multi-speckle measurements [17] by employing a large number of detectors at a given spot on the tissue surface is not feasible for the type of dense sampling [20] that is desired for practical DCT applications.
An alternative would be to employ large arrays of detectors, for example using a photoncounting CCD, to capture the intensity fluctuations of many speckles simultaneously and then computing the intensity autocorrelation. This has been attempted in DWS studies [21,22] but it is not yet feasible for biomedical applications where the dynamics of the speckles, i.e. the decay of the intensity auto-correlation, is much faster than the frame rates achievable by current camera technologies. We note that the speckle contrast measured using a point source illumination with varying coherence was utilized in [23] to measure the optical scattering coefficient rather than blood flow.
In this work, we present a new speckle contrast based tomographic system-speckle contrast optical tomography (SCOT)-that allows us to to probe heterogeneities in the dynamics of turbid media, like blood flow in deep tissues. SCOT is similar to DCS in the sense that it works in the multiple scattering, photon diffusion regime. SCOT merges the deep tissue flow measurement capabilities of DCS and the relatively inexpensive detectors with high frame rates as used in LSF. An improved SNR and a broad field-of-view is achieved by using a 2D detector configuration whereas the simultaneous measurement of large number of speckles (of order a million) over many source-detector separations gives the ability for depth resolution. We describe the theoretical background, derive a perturbation equation based on first Born approximation, describe an inversion algorithm and demonstrate the feasibility of the technique in tissue simulating liquid phantoms.
The rest of the paper is organized as follows. The theory and the inversion algorithm of SCOT are explained in Section 2. We describe the speckle flow imaging with correlation diffusion based photon propagation model and derive the sensitivity relation for the SCOT based on Born approximation. The inversion algorithm for SCOT based on the derived sensitivity relation is also presented. The experimental method of SCOT is presented in section 3 where we describe the optical instrumentation needed to acquire the intensity data for computing the speckle contrast. The reconstruction results are presented in Section 4 and the discussion and conclusions drawn from these studies are presented in Section 5.

Correlation diffusion equation and the laser speckle contrast
The propagation of light in multiple scattering media with static scatterers characterized by the optical absorption (µ a ) and scattering (µ s ) coefficients obeys the radiative transfer equation (RTE) [3]. Under a set of well-defined and validated assumptions , one can arrive at the diffusion equation model [3] for photon propagation which is a a simplified form of RTE. Similarly, in order to study the photon propagation under dynamic scatterers, correlation transport equation (CTE) [24,25] is adopted, which is analogous to the RTE for static scatterers. From CTE, again using the diffusion approximation, correlation diffusion equation (CDE) [10,12,26] is derived which describes the propagation of electric field (E) autocorrelation G 1 (r, τ) =< E * (r,t)E(r,t + τ) > as given in Eq.1: where the reduced scattering coefficient µ ′ s = (1 − g)µ s , g is the anisotropy factor of scattering, D = is the optical diffusion coefficient, and S 0 (r − r 0 ) is the isotropic point source located at r = r 0 . k 0 = 2πn 0 λ is the modulus of propagation vector of light where n 0 is the refractive index of free medium and λ is the wavelength of light. Here r = (x, y, z) is the spatial co-ordinate in three dimensional space and τ is the correlation time.
The dynamics of the scattering medium is modeled by < ∆r 2 (τ) > which is the mean square displacement (MSD) of the scatterers [3]. When the motion of scatterers is modeled as Brownian motion then < ∆r 2 (τ) >= 6D B τ , where D B is the particle diffusion coefficient. Under a directed capillary flow of the scatterers the MSD is given by < ∆r 2 (τ) >= V 2 τ 2 where V 2 is the square of the effective (average) velocity of the scattering particles. We consider here both the Brownian motion as well as the directed flow and hence MSD is given by < ∆r 2 (τ) >= 6D B τ + V 2 τ 2 . Other more general formulations are also possible to account for different types of motion and effects [3,27].
In DCS, the measurement available is the intensity autocorrelation. g 2 (r, τ) =< I(t)I(t + τ) >, instead of the field autocorrelation which is related to normalized field autocorrelation g 1 (r, τ) = G 1 (r,τ) G 1 (r,0) through the Siegert relation [28]: g 2 (τ) = 1 + β |g 1 (τ)| 2 . Here β is a constant, typically less than one, whose value depends on the detection optics as well as the different optical modes in the measurement. Together with the CDE in Eq.1 and the Siegert relation, DCS quantitatively measures the MSD. For diffuse correlation tomography (DCT), in an analogous fashion to DOT, a perturbation approach and an inverse algorithm is utilized connecting the local perturbations of the MSD and the measurement, g 2 (r, τ) for several source-detector pairs (r) and correlation delays (τ). This is then used to recover the spatial distribution of D B or V 2 [11,13,16]. We note that the reported quantity in DCS/DCT and hence in SCOT is an effective diffusion coefficient (D B ) obtained from fitting to correlation data and it is not the traditional thermal Brownian diffusion coefficient of red blood cells. It is, however, well validated that this model represents the absolute and relative blood flow in tissues [3,4,29]. We also note that the use of the reflection geometry is challenging but is feasible and useful as it is the case for any diffuse optical tomography problem including DCT [11,16,30,31] application.
In speckle contrast based measurements, instead of g 2 (τ), another statistical quantity called the speckle contrast, κ, is used for flow measurement. The speckle contrast is defined as the ratio of the standard deviation (σ I ) of measured intensity to its mean (µ I ) value in the spatial domain, i.e. [32], where T is the exposure time of the detection system. κ varies between zero and one and higher values indicate slower fluctuations of the scatterers. In LSF, single scattering approximations and uniform illumination of the sample of interest is utilized relating κ to blood flow at superficial layers. By using the Siegert relation the speckle contrast can be expressed in terms of the normalized field autocorrelation g 1 (r, τ) as [33], We note that this equation with the appropriate Green's function could be used to fit for the background blood flow by using data measured at different source detector separations. The relation connecting κ and g 1 (r, τ) as given in Eq.3 along with the CDE in Eq.1 can be made if we consider point sources and utilize solutions of Eq.1 to obtain g 1 . We now take this approach which generalize the previous works and can be used to generate quantitative measurements. We use the Green's function solution of CDE in Eq.1 for a semi-infinite geometry [16] and the expression in Eq.3 at different source-detector separations and/or exposure times. The use of an analytic Green's function is adopted here for simplicity but more complex methods such as the finite element method (FEM) can also be employed to incorporate complex geometries with inhomogeneous distribution of tissue properties. This was previously reported for both DOT and DCT [13,30] and can directly be applied to SCOT.
The typical behavior of κ with respect to exposure time and the source-detector separation, r, are shown in Figs. 1(a) and (b) respectively. Here we have used µ a = 0.03cm −1 , µ ′ s = 6.31cm −1 , D B = 10 −8 cm 2 /s, β = 0.5 and λ = 7.85 × 10 −5 cm for the computation of Green's function. As expected, the speckle contrast decreases as exposure time is increased for a given sourcedetector separation and also decreases with increasing source-detector separation. It is this dependence on these parameters that we will utilize in SCOT to obtain three-dimensional (3D) images of the distribution of the dynamics of the probed tissue volume.

Born approximation for speckle contrast optical tomography
In order to be able to carry out tomography, we should first develop the inverse imaging problem where the 3D distribution of MSD, ∆r 2 (τ), inside the volume of interest is recovered from the measurement of the two dimensional distribution of intensity speckle contrast, κ, at the surface using the CDE in Eq.1.
We derive the linear Born approximation [19,30,34] to reconstruct ∆r 2 (τ) from κ as explained below. In the future, this could readily be extended to other perturbation methods and non-linear approaches.
The first Born approximation states that scattered field G sc is negligibly small compared to the background (homogeneous) field G 0 1 , i.e., G sc << G 0 1 , and hence the above nonlinear integral equation reduces to a linear integral equation in G sc given by Since κ depends on the normalized field autocorrelation g 1 (τ), in order to compute the perturbations in measurement due to flow term V b , we divide the expression which gives When τ = 0 , G sc = 0 and hence G 0 Hence, the measurement κ corresponding to the inhomogeneous normalized field autocorrelation g 1 under perturbations due to flow is given by In order to simplify the above expression, the condition for first Born approximation (G sc << G 0 1 ) is again invoked resulting in Therefore, the expression for κ reduces to By defining the baseline speckle contrast as κ 2 0 ≡ 2β and ∆κ 2 ≡ κ 2 − κ 2 0 , we proceed to write the change in speckle contrast due the perturbation in MSD ∆r 2 (τ) δ as Here κ 0 is the baseline speckle contrast measured without any perturbation in MSD i.e., when ∆r 2 (τ) δ = 0. We consider a special case of the above sensitivity relation with D δ B = 0 and (V 2 ) 0 = 0 i.e., there is no perturbation in Brownian motion and the perturbation to the system is introduced as a deterministic flow represented by V 2 . We note that the impression of the blood flow in living tissues on the measured quantity is better modelled as ∆r 2 (τ) = D B τ. Here we have adopted a flow imaging experiment which uses ∆r 2 (τ) = V 2 τ 2 . The sensitivity relation in Eq.11 can be employed for living tissue by setting V = 0 and (V 2 ) δ = 0 and considering D 0 B as the baseline flow with D δ B as the perturbation in flow from the baseline value. The estimation of errors associated with the Born's approximation in the context of DOT is addressed in [35] which can be used to estimate the errors for the results presented here, but it is beyond the scope of this paper. .
In order to demonstrate that the so-called banana path of photon propagation is still preserved for the derived Jacobian, in the special case that is considered, we compute the right hand side of the sensitivity relation in Eq.11 in a three dimensional discretized grid (please see Section 2.3 and Fig. 2(b) ). Here analytic Green's function for semi-infinite geometry with homogeneous distribution of tissue properties is used to compute the Jacobian. But the expression for the Jacobian as given in Eq.11 is derived without any assumption on the form of the solution and hence this formalism permits the use of other forward solvers for CDE in Eq.1.
The plot of logarithm of Jacobian in the YZ plane for X co-ordinate=1.8 cm for a source and the detector positioned at (x,y,z)=(1.8 0 2.0) and (x,y,z)=(1.8 1.5 1.08) respectively is shown in Fig. 2(a).

Inversion algorithm
The computation of the sensitivity relation (Eq.11) is done in a discretized slab geometry which corresponds to the size of the object to be imaged in the experiment as well as the scanning pattern of the source relative to detector positions. The computational domain is a rectangular slab of size N x × N y × N z as shown in Fig. 2(b). The distribution of the flow inside the tube is depicted in the three dimensional slice plot in Fig. 2(c) where the geometry shown is used for the simulations. As shown in Fig. 2(d), the sources are scanned along the XZ plane (Y = 0) and the intensity images are collected from the X-Z plane at Y = N y (plane ABCD) which serves as the detector plane.
The first part of the computational algorithm involves the computation of the speckle contrast from experimentally measured raw intensity images. In order to calculate speckle contrast, X-Y co-ordinates of sources from the intensity images were computed using the centre of mass of  images and with the help of a Canny edge detecting algorithm. For each N s sources we have used N d detectors in the detector plane ABCD.
We discretize the (X,Y,Z) co-ordinates into v x , v y and v z points respectively, which gives a total of v x × v y × v z voxels for the three dimensional slab geometry. The baseline speckle contrast (κ 0 ) corresponding to the D 0 B , is computed using the semi-infinite Green's function solution of CDE in Eq.1. The term in left hand side of the Eq.11, which is the perturbation in speckle contrast due to the flow relative to its baseline value (E = ∆κ 2 ), is computed by subtracting the baseline speckle contrast from the speckle contrast computed from intensity images measured in the presence of flow. The right hand side term in Eq.11 is computed using the rectangular geometry with N V = v x × v y × v z voxels and for N s × N d source-detector pairs which gives the Jacobian matrix (J) of size (N s × N d ) × N V .
In particular, we discretize the (X,Y,Z) co-ordinates into 14, 8 and 16 points respectively which gives a total of 1792 voxels with a volume of 0.008 cm 3 for the three dimensional slab geometry. Finally, the Jacobian is normalized to getJ = J • B where ′ • ′ denotes Hadamard product defined as (J • B) i, j = (J) i, j (B) i, j (point wise multiplication). Here B is a matrix of size (N s × N d ) × N V whose rows are the vector b such that b i = 1 a i i = 1...N v and a = (c T + λ 2 max(c T )) [16,36]. Here indices in suffix position is used to denote the elements of the matrix, the vector c, of size N V × 1, is the sum of rows of the Jacobian and 'T ' notates the transpose of the matrix. In the discretized domain, Eq.11 can be re-casted in terms of the normalized JacobianJ which gives where, due to the ill-posedness of the system of equation, we have used λ = λ 1 max(diag(S)).
Here S is the diagonal matrix obtained using the singular value decomposition of the matrixJ. We have found, by trial and error, that λ 2 = 10λ 1 where we have taken λ 1 = 0.1 gives optimal results. The flow velocity is computed by solving the linear system of regularized discrete sensitivity equations given in Eq.12.

Experimental method
The experimental set up is shown in Fig. 3. A temperature controlled continuous laser diode (Thorlabs L785P090, 785 nm, 90 mW) was focused down to a beam of <1 mm diameter to probe the sample. The sample, which consists of a 1% Lipofundin MCT/LCT solution (B.BRAUN, Germany) in water with µ a =0.03 cm −1 , µ ′ s =6.31 cm −1 (both at 785nm) was filled in a transparent plastic container of size N x = 3.8 cm, N y = 1.5 cm and N z = 5 cm as shown in Fig. 2(b). The light source was focused on the bottom of the sample and the produced speckle patterns were imaged from the top with a a monochrome scientific complementary metal-oxide-semiconductor camera (sCMOS; Orca flash4.0, Hamamatsu, Japan) at the top surface of the sample. The horizontal field of view was ≈ 4 cm, resulting in a pixel diameter of 3 × 10 −4 cm. A f/# of 16 was chosen to roughly match the speckle size to pixel size. The exposure time of the camera was set to 1 ms. Each pixel over the image corresponds to a specific distance from the source, and hence, the use of a camera provides dense spatial sampling with a large field-of-view.
A pair of galvo controlled mirrors were used to scan the laser point source in the bottom X-Z plane of the container. We have used N s = 75 different source positions arranged in an array of 3 rows with each having 25 sources. This array is homogeneously distributed in the field of view of the camera (3.9 cm x 3.7 cm in XZ plane ) and each source corresponds to one position of the focused laser beam. The laser was set in every position during 0.5 seconds to acquire 35 intensity images per source, with a 1 ms exposure time. The scanning positions of the source and the detectors are depicted in Fig. 2(d).
As a perturbation, a slightly turbid plastic tube of 0.4 cm diameter is introduced inside the rectangular container through which the same liquid as the rest of the phantom is pumped using a peristaltic pump. Therefore, the optical properties within the tube were same as the background and the tube walls were found to have a negligible effect in the background measured speckle contrast. The following velocities were utilized: V = 0.11, 0.21, 0.32, 0.43, 0.64, 0.85, 1.06, 2.12 and 3.18 (in cm/s). Here the pump gives the flow ( F in cm 3 /s) and we related it to the velocity by using the diameter (d) of the tube as V = 4F πd 2 . For each velocity, we have recorded thirty five transmitted intensity images. The rationale behind choosing a tube with 0.4 cm is to mimic the spatial distribution of the cerebral hemodynamic response to functional stimulation and induced ischemia in rodents [16,37]. The range of velocities is chosen to cover an extreme range of perturbations in the signals.
For the analysis, for each source, N d = 75 detectors were defined, located at XZ plane for Y=1.5 cm thus comprising a total of N s × N d = 5625 source-detector pairs which serves as the SCOT data. κ was calculated for each detector position using a 5x5 pixel window. These values are averaged over 35 images corresponding to each source and using Eq. 2 the speckle contrast for each detector is computed.

Shot noise correction in speckle contrast
In all experiments, the theoretical speckle contrast as shown in Fig. 1(b) will be different from the speckle contrast computed from experimentally measured intensity images due to noise. One of the detrimental noise sources is the shot noise which obey the Poisson statistics. Hence the speckle contrast due to shot noise can be written as κ s = √ γI γI = 1 √ γI , where γ is the ratio of full well capacity of CCD/sCMOS camera to its analog-to-digital conversion bits . The presence of shot noise will result in the speckle contrast computed from experimentally measured intensity images to increase with respect to spatial variable r, in contrary to theoretical behavior as shown in Fig. 1(b). This is evident from the expression for speckle contrast due to shot noise, κ s , which is inversely proportional to intensity and hence directly proportional to r for a point source illumination. In order to reduce out the effect of the shot noise in the speckle contrast we define a corrected speckle contrast [38] which is κ c = (κ 2 − κ 2 s ) where κ s = 1 √ γI with γ = 0.4578 (defined for our specific camera model. The corrected speckle contrast behaves more closer to the speckle contrast derived from theoretical model and hence we use the corrected speckle contrast for SCOT (Please see Fig. 4(a) and 4(b) in section 4 ). The dynamic range could be further extended by improving the detected SNR and taking into account other significant sources of noise that introduce systematic errors.

Validation of the point source model and data
The speckle contrast due to the Brownian motion represented by κ 0 has to be determined apriori for the SCOT inversion procedure. Hence the speckle contrast measurement computed with transmitted intensity images from the Lipofundin phantom illuminated by a point source is fitted against the numerically computed speckle contrast using CDE. In the fitting algorithm, based on nonlinear least square minimization, we have used the experimentally measured values of optical absorption (µ a = 0.03 cm −1 ) and scattering coefficient (µ ′ s = 6.31 cm −1 ) while the algorithm minimizes for D B . Figure 4(a) shows the speckle contrast, computed from experimental data, with and with out shot noise correction. Figure 4(b) shows the computed speckle contrast (from Equation 3 ) fitted against the measured speckle contrast (shot noise corrected) as a function of source distance separation in centimeters. We would like to mention that the speckle contrast up to a source-distance separation of 2.3 cm is only used for fitting since the systematic deviations form the theory due to uncorrected noise factors increases considerably after this separation . The D B = 1.68 × 10 −8 cm 2 /s we obtained by the fitting algorithm is in agreement with the D B = 1.95 × 10 −8 cm 2 /s that we obtained using DCS measurement.  In Fig. 5 we plot the perturbation in the speckle contrast from its background value along the Z-direction as a function of two different flows. Here v 1 is three times higher than the other (v 2 ). We clearly see the velocity dependent change in the speckle contrast due to the flow in the tube. This is essentially a demonstration that the left hand side of Eq.11 is sensitive to the flow velocity in the tube.

Experimental demonstration of SCOT
As noted above, we consider only those detectors that lie within 2.3 cm of each source position. As sources are illuminated in a bounded rectangular geometry, the number of detectors per each source within a distance of 2.3 cm will vary according to the source positions. Therefore, the size of the Jacobian with the new set of detectors is 5269 × 1792. The normalized system of equation, as given in Equation 12, is solved for (V 2 ), whose square root gives the flow velocity. The reconstructed flow profile, for the highest value of flow (3.18 cm/s), as a three dimensional slice plot is shown in Fig. 6(a). Similar plots for the velocities 1.06 cm/s and 0.32 cm/s are shown in Figs. 6(b) and 6(c) respectively. The tube is clearly visible, albeit with a relatively poor resolution as expected from DOT images. Different velocities provide different amounts of contrast. To quantify the observed changes in velocity, we assign a predetermined volume (matching the original position of the tube) which comprises of the rectangular region formed by (X=0.5 to 3.1 cm, Y=0.65 cm to 0.85 cm and Z=1.7 cm to 1.9 cm ) with a total volume of  Fig. 6(d). The normalization is done by dividing the original and reconstructed flow corresponding to the flow value of 0.85 cm/s. A linear fit, using the data from original velocities ranging from 0.11 cm/s to 1.06 cm/s of the reconstructed flow gives a slope of 0.97. This slope is quite encouraging for this limited range. The underestimation for the larger perturbations, i.e. for larger velocities, is presumably due to the failure of the linearized Born approximation [39].

Discussions and conclusion
We have presented a new tomographic imaging method, speckle contrast optical tomography (SCOT), to measure the blood flow distributed in deeper regions of the tissue. We have proposed and demonstrated that by scanning the point laser source over the sample and acquiring multiple speckle measurements simultaneously using a array of detectors like a CCD or CMOS camera, we are able to model the speckle statistics with a photon diffusion model, fit the model to the data to obtain information about sample dynamics and use in a tomographic inverse problem.
After deriving the physical model and a linearized inversion model, we have carried exper-iments in a liquid phantom with a cylindrical inclusion with varying flow rates of scatters in the transmission geometry. We were able to obtain three dimensional reconstructions and quantify the flow rate over a large range. For the in vivo applications, this technique would allow us to image local or temporal variations in blood flow for example due to neuronal stimuli, pharmacological or physiological changes in time or heterogeneities due to local ischemia or a tumor.
In general, the inverse problem associated with the recovery of MSD from the speckle contrast measurement can be re-casted as nonlinear optimization problem. We have adopted the standard procedure to linearize the non-linear inverse problem using the Born approximation. We have derived the sensitivity relation connecting the perturbations in speckle contrast measurements from its baseline value to the changes in flow velocity. In the discretized geometry, the sensitivity relation is called the Jacobian which was plotted for a given source-detector pair to show that it preserves the so-called banana path as observed in DOT. We have further plotted the corrected speckle contrast as a function of spatial co-ordinates versus different values of flow velocity to show the sensitivity of the data for inversion to the perturbations in flow. Having computed the Jacobian and measure the perturbation in speckle contrast from its baseline value, we presented the recovery of the flow velocity using Tikhonov-regularized least square minimization. In a nutshell, SCOT overcomes the limited dynamic range and low SNR of DWS (and hence DCT) and lower penetration depth of LSF and provides a three dimensional distribution of flow in tissue using faster and cheaper instrumentation.
In past, various proposals were made to improve the depth of imaging. One study [40] proposed a line beam scanning illumination instead of the uniform illumination. However, while promising, this study did not develop a physical model that allowed a quantitative measurement. In another study [41], the authors have used a transmission geometry through the finger joint, 1.0-1.5 cm, but did not use a proper physical model or achieved any depth resolution. Depth sensitive speckle contrast measurement using a sinusoidally modulated source is used in [42], where a pair of confined flows separated 2 mm apart and at a depth of 2 mm and 4 mm were measured. Although this work presents the relation connecting speckle contrast to the electric field autocorrelation for sinusoidal source excitation, a quantitative recovery of the effective flow velocity from the measured speckle contrast was not presented.
Furthermore, we note that quantitative measurements with point-sources of homogeneous media are feasible, i.e., this work also demonstrates a speckle contrast optical spectroscopy (SCOS) that is discussed in details elsewhere [?] which is the basis of the model utilized for the SCOT method. The SCOS method is a more general and quantitative method compared to previous point source speckle contrast studies [43,44]. In a recent study [43], authors measure flow in tissue mimicking phantoms as well as in human arm using speckle contrast measurements but with point sources and detection at a distance. This study was further extended [44] by the same group to multi channel deep tissue flowmetry. However, they did not extract the dynamic parameters from the laser speckle contrast data, but, instead, they have worked in a range where the two are linearly related and considered the inverse of the laser speckle contrast as an index of blood flow. These studies have shown that quantitative laser speckle contrast measurements with point sources may be feasible by introduction of a proper inverse model relating the dynamics of the scatterers to the speckle contrast. Furthermore, these studies have assumed a homogeneous distribution of blood flow in the phantoms or tissues and did not attempt to reconstruct the distributions of the variations in the dynamics.
Previously, deep tissue blood flow measurements was made possible by diffuse correlation spectroscopy (DCS) [3,10,26] which relies on the measurement of the intensity autocorrelation of the random fluctuations of multiply-scattered light intensity. DCS is essentially a unification of the multiple scattering models for LSF and diffuse wave spectroscopy (DWS) in the setting of the tissue optics. A coherent laser point source is used in DCS which enables us to probe several centimetres of tissue with high temporal (≈ 100ms) but limited spatial(≈ 0.1 cm) resolution. Traditional DCS measurement employs individual detectors to capture the temporal statistics of intensity from each speckle [11,12,45].
We acknowledge that the spatial localization of the reconstructed flow velocity values needs further improvement as evident from the Figs. 6(a)-(d). This suggests us to use a denser scanning of the the source in the XZ plane as the slices near to the source positions are having good spatial localization of the flow profile. This also reduces the inherent ill-posedness associated with the inverse problems involving the diffusion equation. This type of denser scanning would have been prohibitively expensive and/or time consuming with the current DCT technologies. However, with SCOT, this can readily be implement.
The limited linear range in measuring the flow velocity as revealed by the plot in Fig. 6(d) suggest to adopt an iterative Born inversion instead of the first Born approximation. As discussed above, the estimation of errors associated with the Born's approximation was previously carried out in the context of DOT [35] and was experimentally demonstrated [39]. A more accurate nonlinear-iterative reconstruction algorithm can be adopted to improve the accuracy of the reconstruction as was already presented for DCT and DOT using numerical methods such as FEM [13,30].
Our results are bit limited in dynamic range since the data is plagued by systematic deviations from the theoretical model due to unaccounted sources of noise. We have considered only shot noise along with a method to remove its effect in the present work whereas other noise sources like read noise, flat pixel noise, variability of hot pixels have to be accounted to derive an accurate noise model. The next step is to translate these measurements to non-invasive, noncontact measurements of the 3D blood flow in vivo.
Deep tissue perfusion imaging is considered as one of the important imaging problems in the medical imaging research. In this work, we have developed and experimentally demonstrated a new optical method for three dimensional imaging of blood flow in tissues. Speckle contrast optical tomography (SCOT) is based on laser speckle contrast approaches and combines the physical models utilized in diffuse correlation tomography (DCT) with laser speckle flowmetry (LSF). The main potential advantage of SCOT is that it combines the deep perfusion imaging capability of DCS along with the rapid data acquisition possible in LSF even with a very dense sampling of sources and detectors without being prohibitively expensive. This makes SCOT a promising technique for deep perfusion imaging.