Quantitative, depth-resolved determination of particle motion using multi-exposure, spatial frequency domain laser speckle imaging

: Laser Speckle Imaging (LSI) is a simple, noninvasive technique for rapid imaging of particle motion in scattering media such as biological tissue. LSI is generally used to derive a qualitative index of relative blood flow due to unknown impact from several variables that affect speckle contrast. These variables may include optical absorption and scattering coefficients, multi-layer dynamics including static, non-ergodic regions, and systematic effects such as laser coherence length. In order to account for these effects and move toward quantitative, depth-resolved LSI, we have developed a method that combines Monte Carlo modeling, multi-exposure speckle imaging (MESI), spatial frequency domain imaging (SFDI), and careful instrument calibration. Monte Carlo models were used to generate total and layer-specific fractional momentum transfer distributions. This information was used to predict speckle contrast as a function of exposure time, spatial frequency, layer thickness, and layer dynamics. To verify with experimental data, controlled phantom experiments with characteristic tissue optical properties were performed using a structured light speckle imaging system. Three main geometries were explored: 1) diffusive dynamic layer beneath a static layer, 2) static layer beneath a diffuse dynamic layer, and 3) directed flow (tube) submerged in a dynamic scattering layer. Data fits were performed using the Monte Carlo model, which accurately reconstructed the type of particle flow (diffusive or directed) in each layer, the layer thickness, and absolute flow speeds to within 15% or better.

An important modeling advancement for DLS in turbid media is the idea that the electricfield temporal autocorrelation function G 1 (τ) = <E(0)E*(τ)> obeys a transport equation [17]. Thus, the autocorrelation function can be found as a function of source-detector separation, for example, using familiar transport solver techniques such as diffusion or Monte Carlo [10]. Making use of the Siegert relation, the field autocorrelation function can be related to the experimental realizable intensity autocorrelation in order to determine the dynamics of light scattering structures [12]. Measurement of the intensity autocorrelation is relatively straightforward for fiber based techniques such as DWS, where the temporal resolution is fine and the distance between the source and detector is clearly defined. In contrast, image based measurements such as LSI generally require longer exposure times and ill defined sourcedetector elements. However, by performing radial Hankel Transforms on results from transport models, one can instead find autocorrelation as a function of spatial frequency [18]. The planar nature of Fourier basis functions are naturally well suited for two-dimensional imaging.
Structured light concepts have been applied to the transport of light intensity as the foundation for a reflectance imaging technique called Spatial Frequency Domain Imaging (SFDI). SFDI has been well established to determine the absorption, scattering, and fluorescent properties in turbid media [19][20][21][22][23][24][25][26]. Using structured light also allows for a tomographic reconstruction of layered, heterogeneous optical properties [22,27], because the depth penetration of the light depends on the spatial frequency of projection.
In a previous manuscript, we applied experimental methods and modeling precedent from SFDI to the transport of correlation by projecting structured coherent light on dynamic scattering objects, and measuring the speckle contrast as a function of spatial frequency [28]. We showed that the speckle contrast follows expected behavior as a function of spatial frequency [29], and that correctly modeling the effect using correlation transport allows one to fit quantitatively for the flow properties of light scattering particles in a homogenous medium [28].
In a biological setting, samples are often layered and heterogeneous. Thus, planar speckle images are formed from photons that may interact with scatterers of different flow dynamics. This issue has been recognized and addressed for samples composed of a static scattering layer above a dynamic layer [30]. Here, the presence of a static top layer introduces nonergodic speckle fluctuations, which affect fitted flow values. To rectify non-ergodic data, the field autocorrelation was re-derived assuming a superposition of fields from photons scattered by static and dynamic particles [15,30]. At least one additional image exposure must be taken to separate the dynamic field component. This method has been used successfully in phantoms and in-vivo [31,32], but the method falls short in some aspects: it does not account for photons that scatter partially in both static and dynamic regions, and it cannot determine the presence of static scatterering layers beneath the sample surface nor situations where both layers are dynamic, but with different flow rates or types (e.g. a large vessel buried beneath a capillary network). Our manuscript demonstrates a method that overcomes these limitations by incorporating multiple exposures and spatial frequencies into the speckle contrast measurement. With this method, we are able to determine quantitative flow properties in multi-layered media.

Theory
LSI measures speckle contrast, which is defined as: where K is the speckle contrast, σ is the intensity standard deviation, and <I> is the mean intensity. This calculation is generally performed using a sliding window filter across the raw speckle image, which then generates a speckle contrast image. For a speckle pattern with some finite integration time T, K is related to the electric field autocorrelation function G 1 (τ) through the expression [14,33] ( ) ( ) Here β represents a multiplicative reduction in contrast not associated with the dynamics of the scatterers, such as depolarization, coherence length, or mismatch between the image pixel and speckle size [34,35]. For turbid media, a correlation diffusion equation can be used to find G 1 analytically as a function of optical properties and the mean square displacement <Δr 2 (τ)> [10]. However, we and others have shown that when imaging speckle over a wide field, the diffusion approximation breaks down more readily because of shorter average photon path lengths and relatively long camera exposure times [28,30]. Therefore, we utilize Monte Carlo (MC) simulation techniques to determine the autocorrelation function. As proposed in [9,10,36], the autocorrelation can be found as a function of separation from a point light source, ρ, by tracking the momentum transfer distribution at discrete ρ bins and numerically integrating over the single scattering decorrelation exponential [37]. The minimum bin resolution determines the maximum spatial frequency able to be modeled without aliasing. In our case we use Δρ = 0.1 mm, which according to Nyquist yields f max = 1/2Δρ = 5 mm −1 , easily sampling our intended maximum spatial frequency of 0.25 mm −1 : Here P(Y) is the distribution of dimensionless momentum transfer, Y = 1-cos(θ), k 0 is the wave number inside the medium, and <Δr 2 (τ)> is the mean-squared displacement. Because MESI affords us the opportunity to fit for the type of motion (diffusive or directed flow), we choose to write the displacement in terms of its flow coefficient and temporal power: Handling multiple layers adds extra difficulty because the decorrelation exponential is a nonlinear function of both the total momentum transfer Y and also the fractional momentum transfer distribution in each layer for that Y. The autocorrelation becomes a multi-integral, which for a two layer (or three layer with identical 1st and 3rd layer) system can be written using conditional probabilities: Here y is the fraction of Y originating from the first layer, d is the first layer depth, and V 1 , V 2 are the flow coefficients of the first (or first and third) and second layer, respectively. Note that the conditional probability distribution P(y|ρ) must be generated for each value of the source-detector separation ρ i (or momentum transfer). We choose a two-layer system as the simplest case, but the method can be extended for arbitrary numbers of layers. However, each additional layer increases the number of conditional probability distributions required by another power.
It is important to realize that the approximation for speckle contrast in the presence of static scatterers, as presented in [10,30,32,34], can be re-derived using our formalism by assuming that the conditional probability is separable and equivalent to P(y|ρ) = wδ(0) + (1w)δ (1), where the weighting factor w is the usual ratio of dynamic to static intensity <I f >/(<I f > + <I s >). This is equivalent to assuming that the speckle contrast arises from photons that scatter only in static or dynamic layers, but not both. Clearly, numerous photons scatter multiply through both layers, hence we feel the full distribution must be calculated and integrated to increase accuracy of the model.
To find the autocorrelation as a function of spatial frequency, one can apply a numerical Hankel Transform to G 1 [28] ( ) ( ) ( ) Finally, the speckle contrast modulation transfer function K(f), can be found by inserting G 1 into Eq. (2). K(f) is the expected speckle contrast for each spatial frequency when the coherent source has spatial structure. Note the order of Eq. (4)-(5) is arbitrary, as the Hankel Transform can be applied to P(Y,ρ), then integrated. Note also the derivation here is done in terms of continuous variables and integrals, but were performed numerically using results provided by Monte Carlo methods. Potential parameters that can be found when fitting to combined multi-exposure and multi-frequency data include: V 1 , V 2 , d, n 1 , n 2 .

Methods
A multi-layer C# based Monte Carlo application was downloaded from the Virtual Photonics Technology Initiative (http://www.virtualphotonics.org/). This application tracks photon weight and momentum transfer as a function of layer and optical properties using standardized Monte Carlo RTE solver methodology [38]. Histograms were generated using one million photons, at appropriate optical properties, for the total and fractional momentum transfer, P(Y,ρ) and P(y|ρ), respectively. These results were Hankel Transformed and numerically integrated using Matlab (Mathworks, Inc) software to give predicted K(f).
Validation experiments were performed using a structured light imaging system shown in through a series of focusing optics onto the sample surface. The optics kept the field of projection constant, and the frequency was changed by translating the stage along the divergent source axis. Images were captured with a Retiga EXI cooled 12-bit CCD camera (Qimaging, Inc), equipped with a 200 mm focal length lens set at F16, giving a field of view approximately 8 mm by 6 mm. The optical imaging axis was angled slightly (~10 degrees) to avoid specular reflection. A cross polarizer was also placed in the camera lens to help further reject specular light and increase speckle contrast.
To determine speckle contrast as a function of spatial frequency we take advantage of the demodulation scheme described in [20]. This technique allows one to isolate the intensity I (or standard deviation) amplitude of a single frequency component by combining images at three equally separated phases φ 1,2,3 = (0, 2π/3, 4π/3) such that: The amplitude of the standard deviation maps can be calculated in the same manner as it scales linearly with intensity. The contrast is then simply K AC = σ AC /I AC . Standard deviation and mean maps were calculated using a small sliding window filter of 7x7 pixels. Border areas within 7 pixels of the edge were discarded. A window size of 7x7 is typical in the literature, and chosen as a reasonable compromise between spatial resolution and increased sampling statistics [39]. All measurements were carefully designed to meet or exceed Nyquist criteria and minimize sampling errors [35]. The value of β was found empirically using a static diffusive target and calculated to be β = 0.81. A variety of static and dynamic phantoms were created to validate our models. Static phantoms were created at varying thicknesses using silicone PDMS. Titanium dioxide was used as a scattering agent, and India ink provided absorption. Phantom thickness was set by shallow molds and verified using calipers. The optical coefficients for static phantoms at 685 nm light were set at approximately μ s ' = 1.2 mm −1 , μ a = 0.005 mm −1 , based on known values for these substances [40][41][42]. Dynamic liquid phantoms were created using polystyrene microspheres of diameter 800 nm (Spherotech, Inc.) as a scattering agent, and India ink. Using the Stokes-Einstein equation for diffusion of spheres in liquid, the diffusion coefficient was calculated to be D b = 6.1x10 −7 mm 2 /s. The scattering coefficient for liquid phantoms at 685 nm was set at μ s ' = 0.6 mm −1 , using an appropriate concentration of microbeads. This concentration was determined by theoretical Mie scattering calculations, conveniently provided by the Oregon Medical Laser Center (http://omlc.ogi.edu/calc/mie_calc.html). Absorption was set at μ a = 0.005 mm −1 .
Phantom experiments were conducted in three separate test configurations shown in Fig.  1.  First, thin static phantoms were placed over a semi-infinite well of dynamic phantom. The layer thickness was varied from 0.4 to 4 mm. This setup was designed to be comparable with previous experiments using MESI to account for static speckle layers [30].
Second, a thin liquid phantom was poured into a container above a semi-infinite layer of static phantom. The layer thickness was varied from 0.75 mm-3.75 mm. For thin liquid layers under 1 mm, a small amount of dish soap (~1 drop) was added to reduce surface tension and decrease beading. The soap was found to have insignificant impact on solution viscosity and optical coefficients. This particular setup was designed to illustrate where the MESI method falls short due to its reliance on static speckles in each image remaining constant with exposure time. For a dynamic top layer, nearly all photons interact with dynamic scatterers but are still significantly influenced by static sections below.
Third, a dynamic scattering liquid was confined in a clear polyurethane tube of radius 2 mm, pumped at a known flow rate, and submerged in a container full of an identical scattering liquid (at rest). The flow rate was controlled with a syringe pump (New Era Pump Systems Inc.) and varied from 0.375 mm/s-1.88 mm/s. The depth of the submerged tube was fixed at 1.25mm. MESI methods designed to account for static speckle contrast would also fail in this situation since both layers are dynamic. The (relatively) slow Brownian diffusion of the top layer dilutes the overall speckle contrast considerably when compared to the fast, directional flow in the tube.
Experimental data for all configurations were collected at 10 exposure times evenly spaced between 0.5 and 14 ms. These values were chosen to maximize the "S" sensitivity curve of speckle contrast to exposure time at our dynamics of interest (i.e. Brownian Motion and Flow < 2 mm/s), as predicted by literature values and forward simulation [30,32]. For configuration one, six spatial frequencies were imaged between 0 and 0.25 mm −1 . For configuration two and three, five spatial frequencies were recorded between 0 and 0.12 mm −1 . In each case these frequency bounds were chosen to maximize sensitivity of the speckle contrast measurement, with higher frequencies eliminated as they approach an asymptote in the momentum transfer distribution. The number of frequencies and exposures was chosen as a compromise between time of measurement and accuracy. Note the minimum number must exceed the number of unknowns, while the optimal number is still unknown and requires further study. After images were recorded, speckle contrast data was averaged across the image at all frequencies for phantom configurations one and two, and then fit as a single data point. Fits were performed using a non-linear least squares solver for combined multi-exposure and spatial frequency data. Four parameters were fit: velocity coefficients V 1 , V 2 , the temporal displacement powers n 1 , n 2 , and the top layer thickness d. For configuration 3, speckle contrast data was binned to a size of 39 by 29 pixels and fits performed for each individual pixel. Sub-regions of interest were then chosen and the fits averaged. Fits were performed in two steps: first, multi-exposure high spatial frequency content (>0.1 mm −1 ) was used to find V 1 and n 1 outside the tube region, where homogeneity can be assumed. Here the shorter mean path lengths of high frequency content served to reduce partial volume effects caused by the tube. Once V 1 was fit, this variable was fixed in the least squares solver and planar multiexposure content (0 mm −1 ) was used to fit for V 2 , n 2 , and d.

Results
Momentum transfer probability distributions as a function of spatial frequency, generated by our Monte Carlo model for μ s ' = 0.6 mm −1 , μ a = 0.005 mm −1 , are shown in Fig. 2. The probability distributions illustrate significant change in momentum transfer and depth penetration with spatial frequency, which is the driving force behind our technique. This dependence is used to fit speckle contrast at multiple exposures and spatial frequencies to quantitative, layered flow values. Figure 3 shows experimental speckle contrast with corresponding fits for phantom configurations one and two. Fig. 3. Validation phantom speckle contrast fits for configurations 1 and 2. Experimental speckle contrast measurements (circles) and subsequent fits (solid line) for five values of thickness d, linearly spaced between 0.7mm (lightest) and 3.0mm (darkest). Notice the drastic difference between the shape of the speckle contrast curve, as low frequencies tend to average both layers, and high frequencies tend toward the dynamics of the top layer. Figure 4 shows fit values corresponding to phantom setups one and two. Note that temporal displacement power n was fit as n = 1 for both cases. Thus, the flow coefficient here is equivalent to six times the Brownian Diffusion Coefficient D b . While the static layer diffusion coefficient is expected to be zero, this value cannot be reached by the fitting algorithm once it becomes arbitrarily small and ceases to affect the contrast. We define this "effective zero" as the point when there is < 3% difference between the effective and actual zero speckle contrast. The effective zero depends on layer thickness.  Recovered n values were found to be one in both cases, and thus we define the flow coefficient as six times the Brownian Diffusion Coefficient D b . Multi-layer fits for dynamic phantom D b (circles) and static phantom (squares) are shown alongside blind single layer fits (x's) and expected values (solid and dashed lines). Blind fits can be seen to vary significantly, over two orders of magnitude in some cases, from expected values. Notice the effective zero line (dashed) depends on layer thickness and depends heavily on the phantom geometry. Depth fits generally underestimate true values by a small margin.
An exploration into partial volume effects, and subsequent reasoning behind the two-step fitting process for phantom configuration three is shown in Fig. 5. Here the speckle contrast image at f = 0 mm −1 is shown, and Region of Interest (ROI) highlighted in a portion of the image expected to be homogenous. Mean speckle contrast in this ROI is calculated for f = 0 mm −1 and f = 0.12 mm −1 , and compared with measurements made on a truly homogenous sample with no tube present. As can be seen from Fig. 5, at f = 0 mm −1 speckle contrast differs substantially from a true homogenous sample, and that difference increases with tube velocity. This is due to the fact that at low projection frequencies, the relatively long average optical path lengths allow a significant fraction of the detected light to interact with the tube [20]. In contrast, high projection frequencies reduce the mean optical path length, and there is essentially no detectable contribution from the deeper dynamic structure. Thus, to increase accuracy and effective resolution, flow coefficients and flow types for the surrounding medium, V 1 and n 1 , are fit using exclusively high frequency data.  . Experimental speckle contrast image for DC spatial frequency (0 mm −1 ), and representative ROI for "homogenous" regions (left). The region appears uniform but is actually highly dependent on tube velocity for low spatial frequencies, as shown in the plot (right). Here, the DC frequency (x's) varies by as much as 13% from a truly homogenous sample. The AC (circles) spatial frequency (0.12 mm −1 ), in contrast, is not heavily influenced by tube velocity due to lower average path lengths. Figure 6 shows multi-exposure fitting results at a range of tube velocity values. Fitted speckle contrast for both "homogenous" sections and high flow "tube" sections of each image are shown in Fig. 6(b). Fit variables V 2 , n 2 , and d, are shown in Fig. 6(c), 6(d), and 6(e), and performed using multi-exposure data at f = 0 mm −1 and previous high frequency fits for V 1 and n 1 , described above. The fits are averaged for the two given ROI's.  Note that in more complicated samples, there may not be a clearly identifiable "homogenous" region, but the operating principle is still exceedingly useful. Essentially, high spatial frequencies gate diffuse photon components from the measurement, and allow for the isolation of surface flow. We believe this information to be valuable in any speckle image fitting process.

Discussion and conclusion
A method for fitting speckle contrast to quantitative flow in layered geometries was shown, which combined features of Monte Carlo modeling, multi-exposure speckle imaging, and spatial frequency domain analysis. To model this process correctly, a key data feature must be tracked in our Monte Carlo algorithm: layer-specific fractional momentum transfer distribution. This distribution shows a dependence of "momentum transfer penetration depth" on spatial frequency, where higher spatial frequencies penetrate less into a scattering medium. Useful parameters such as flow coefficients and temporal displacement power also depend strongly on exposure time. Thus, we structured coherent light and measured contrast at several exposures to take advantage of these complementary effects. Using controlled phantom experiments with known flow values in three different geometries, we showed that all parameters of interest could be fit accurately to speckle contrast maps.
Translating this method to biological samples will admittedly present extra complications, as these are typically composed of multiple intersecting layers with highly variable flow speeds. However, incorporating multiple spatial frequencies and exposure times into speckle imaging provides researchers necessary additional elements of control in the sampled photon path length and sensitivity to flow speeds present in a sample volume. By investigating the potential fitting power behind this additional information content, our manuscript represents a building block towards biological applicability. We believe that with the advent of more advanced imaging sensors, lasers, and Monte Carlo solvers, researchers will be able to quantitatively image flow in a substantially broader range of conditions.
In ongoing research, we plan to expand our technique to in-vivo tissue studies. Important applications include measuring blood flow and metabolism in pre-clinical rodent brain models, retina, and skin. In each of these cases the presence of layered structures with varying optical properties can have a sizeable impact on the measured speckle contrast. We hope to combine optical absorption tomography and layered speckle measurements for a more accurate rendering of metabolic tissue function.