Quantitative determination of dynamical properties using coherent spatial frequency domain imaging.

Laser speckle imaging (LSI) is a fast, noninvasive method to obtain relative particle dynamics in highly light scattering media, such as biological tissue. To make quantitative measurements, we combine LSI with spatial frequency domain imaging, a technique where samples are illuminated with sinusoidal intensity patterns of light that control the characteristic path lengths of photons in the sample. We use both diffusion and radiative transport to predict the speckle contrast of coherent light remitted from turbid media. We validate our technique by mea- suring known Brownian diffusion coefficients ( D b ) of scattering liquid phantoms. Monte Carlo (MC) simulations of radiative transport were found to provide the most accurate contrast predictions. For polystyrene microspheres of radius 800 nm in water, the expected and fit D b using radiative transport were 6 : 10 E – 07 and 7 : 10 E – 07 mm 2 = s , re-spectively. For polystyrene microspheres of radius 1026 nm in water, the expected and fit D b were 4 : 7 E – 07 and 5 : 35 mm 2 = s , respectively. For scattering particles in water – glycerin solutions, the fit fractional changes in D b with changes in viscosity were all found to be within 3% of the expected value. © 2011 Optical Society of America OCIS codes: 110.6150, 170.6480.


INTRODUCTION
Dynamic light scattering (DLS) is a method that utilizes coherent light to estimate the motion of particles in scattering media. DLS methodology is used across scientific disciplines to characterize samples in physical, chemical, and biological settings [1][2][3]. In biomedical applications, several DLS methods are used, including: laser Doppler flowmetry [4], diffusing wave spectroscopy (DWS) [5], and laser speckle imaging (LSI) [6][7][8]. In all of these, the decay of the autocorrelation function of coherent light is related to the degree of scatterer motion.
LSI is a particularly attractive method because of its ability to image wide fields of view quickly in a simple and inexpensive manner. However, quantitative measurements using the backscattering geometry of LSI have remained elusive because the decay of the autocorrelation function depends on photon path length and is thus dependent on the absorption and scattering coefficients of the medium. DWS overcomes this difficulty via use of fixed source-detector separations and diffusion-based models of photon propagation [9]. However, due to the limited number of detectors used with DWS, high-resolution images of blood flow are not possible, in contrast to the megapixel spatial sampling commonly used with CCD-based LSI. In addition, LSI is used to image tissue near boundaries and measure blood flow in highly absorbing blood vessels, conditions for which the diffusion approximation breaks down. In order to overcome these limitations, we make the following modifications to LSI. First, we model LSI using MC simulations of the correlation transport equation. Second, we introduce a technique, coherent spatial frequency domain imaging (cSFDI), that allows us to control the characteristic path lengths of photons in our LSI experiments.
Conventional spatial frequency domain imaging (SFDI) uses incoherent light. SFDI is a fast, noncontact optical technique that projects sinusoidal patterns of light on a sample surface. By modeling the propagation of the patterned light through the medium and measuring the backscattered reflectance at specific spatial frequencies, one can image the absorption and scattering properties of the sample using a CCD camera. The technique is described in detail in [9], and has since expanded its scope in terms of application and theory [9][10][11][12][13][14][15]. Briefly, the low-pass filter characteristics of a turbid medium are determined by measuring the modulated transfer function (MTF) of the diffusing light. The attenuation is isolated on a pixel-by-pixel basis as a function of spatial frequency using a multiphase demodulation scheme, discussed in Subsection 3.B. Once the MTF is found, it can be fit to a light-propagation model to determine a unique pair of absorption and scattering coefficients.
In this paper we demonstrate that, by modulating coherent light in a similar fashion, one may also measure speckle contrast as a function of spatial frequency. We demonstrate a distinct advantage to this over standard LSI because, by changing the spatial frequency of the projected patterns, we are able to control the characteristic photon path lengths, and thus the speckle contrast. We can then fit the measured speckle contrast, as a function of projection frequency, to a light-propagation model and quantitatively determine the rate of motion for scattering particles. The approach is validated with experiments on liquid phantoms of polystyrene spheres and fatty emulsions in glycerin-water solutions. Light propagation is modeled with both analytic solutions to the standard diffusion approximation (SDA) and MC simulations of radiative transport. We show that MC modeling of radiative transport gives much more accurate predictions of speckle contrast than the diffusion approximation.

METHODS
A. Laser Speckle Imaging For an ergodic system, one measure of the correlation manifests itself in the speckle contrast of an image. Speckle refers to the granular intensity pattern, caused by interference, that appears when coherent light illuminates a scattering object. The contrast of this pattern is defined as where K is the contrast,σ is the intensity standard deviation, and hIi is the mean. In the case of static scatterers or rough reflecting surfaces, the probability distribution of the intensity follows a negative exponential shape, and the contrast becomes unity [16]. If the scattering objects are in motion, the contrast will be reduced by a quantity dependent on the rate of motion and magnitude of momentum transfer at each photon collision. In LSI, this contrast analysis is performed in small sampling windows several pixels wide to preserve spatial resolution of blood vessels and other structures. Goodman, Bandyopadhyay, and others have derived the contrast of integrated speckle intensity given a uniform integration window as [16,17] where K is the contrast, T is the integration time, and g 1 is the normalized electric field autocorrelation function defined as g 1 ðτÞ ¼ hEð0ÞE Ã ðτÞi=hjEð0Þj 2 i. Here the Siegert relation has been employed to transition from field statistics to intensity [18]. Maret and Wolf then derived the electric field autocorrelation of a single photon as it travels through a scattering medium [19]: where each i is a scattering event, k 0 is the wavenumber inside the medium, k 0 ¼ 2πn=λ, hΔr 2 i is the mean square displacement, and q is the momentum transfer jqj ¼ 2k 0 sinðθ=2Þ. The particle rate of motion information is contained in hΔr 2 i, and depends on the nature of the dynamics in question. For diffusive dynamics, such as Brownian motion, where D b is the diffusion coefficient. For directed motion, hΔr 2 i ¼ ðV τÞ 2 where V is the velocity. In LSI, hΔr 2 i is fit or solved for, and typically the dynamics are assumed a priori. Because the pixel intensity will be a summed contribution of multiple photons, the amount of momentum transfer follows a probability distribution dependent on the optical properties. For such distributions, there are several ways to derive G 1 depending on the manner in which it is assumed light propagates through the medium. The autocorrelation function, which depends on the dynamics of the scattering particles, can then be related to the measured quantity, speckle contrast, using Eq. (2). One method for obtaining the autocorrelation function G 1 is to assume that speckle is a single scattering phenomenon at known angle of incidence, in which case the value of jqj can be found, reducing Eq. (3) to a simple exponential. For relative changes in flow, fitting to this expression often suffices, and has the advantage of being simple enough to be performed in real time [8]. However, it is well known that speckle patterns produced from a turbid medium is a multiple scattering phenomena [20]. To achieve quantification, we instead focus on a model for the correlation that assumes multiple scattering.
In DWS, a sufficiently large number of scattering events is assumed such that the sine of each scattering angle can be replaced by the average. The total momentum transferred is then just the total number of scattering events multiplied by this average. For a single photon, this reduces the autocorrelation function to where s and l Ã are the total path length and mean free path of photon propagation, respectively, and their ratio gives the total number of scattering events. The breakdown of this assumption for short path length photons has been explored in some detail [21]. When extending this result to many photons, one now only needs to integrate over the distribution of path lengths PðsÞ. Using the standard diffusion approximation, this distribution has been calculated as a function of source-detector separation jr − r s j and the integral can be done analytically [22]: It is important to note that this result has also been derived in the context of correlation transport. If correlation is the transporting quantity, one can use the SDA to reduce the equation to [20] where SðrÞ is the source, μ tr ¼ μ a þ μ 0 s , μ eff ðτÞ ¼ ð3μ a;dyn ðτÞμ 0 tr Þ 1=2 , and μ a;dyn ðτÞ ¼ μ a þ 1=3⋅μ 0 s k 2 0 hΔr 2 ðτÞi. This is identical to the diffusion of the photon density, but with the absorption coefficient replaced by this new "dynamic" absorption term μ a;dyn . If a point source is assumed, this can be solved to obtain the correlation in Eq. (5). The dynamic absorption is useful because it allows one to draw on precedents with diffusing flux, but it also adds new restrictions on the validity of the approximation. These issues will be discussed below.

B. Correlation Diffusion in the Spatial Frequency Domain
While obtaining flux and correlation as a function of separation from a pointlike light source is advantageous in optical fiber methods such as DWS, LSI is performed with planar illumination where such a source does not exist. However, recently, imaging techniques have been developed for sources that are two-dimensional plane waves at varying spatial frequencies [9,12,14,23,24]. This can also be viewed as imaging in the Fourier inverse domain of spatial separation, spatial frequency, where the projection patterns are the Fourier basis functions. Analogous to the way that DWS takes advantage of correlation at multiple source-detector separations, we use discrete spatial frequencies of projected light to increase information content and help achieve true quantitative measurements. In diffusion-based SFDI, the source term in the diffusion equation [Eq. (6)] is replaced by a function of spatial frequency. For a source SðkÞ ¼ sinðkx þ φÞ, the solution for the correlation of remitted light is also sinusoidal, with amplitude given by [9] Here the mean square displacement appears inside the effective absorption term, and can thus be fit using this analytic expression.

C. Correlation Monte Carlo in the Spatial Frequency Domain
As opposed to the analytic diffusion model, photon correlation can also be found numerically with MC techniques. Using well-developed MC algorithms for photon propagation through a turbid medium [25,26], one can record the momentum transfer of each photon-scatterer collision directly, and create a numerical distribution. This avoids the DWS assumption that total momentum transfer can be replaced by the number of scattering events. The autocorrelation can then be calculated as [21,27,28] where PðY Þ is the distribution of dimensionless momentum transfer Y ¼ 1 − cosðθÞ. This integral must be evaluated numerically and can be viewed as a weighted average of the exponential function.
By taking a Fourier transform of autocorrelation function determined with Eq. (8), we find the correlation as a function of spatial frequency. Since the problem is spherically symmetric, we calculate the spatial frequency components with a Hankel transform: A simple time-resolved MC photon propagation algorithm was obtained from the Oregon Medical Laser Center website [29]. An additional feature was added to histogram the weighted, dimensionless momentum transfer. One million photons were propagated through a semi-infinite medium to obtain PðY Þ at 100 discrete radial bins, linearly spaced between 0 and 20 mm. The autocorrelation function was then determined by using Eqs. (8) and (9).

EXPERIMENT A. Instrumentation
The experimental setup is illustrated in Fig. 1. Light from a coherent 30 mW He-Ne laser source (Edmund Optics) was expanded and reflected off a liquid-crystal-on-silicon (LCOS) (Holoeye Inc) spatial light modulator (SLM) and directed onto the sample. The remitted speckle pattern was recorded using a 12 bit thermoelectrically cooled CCD camera (Qimaging Inc). The SLM was programmed to project light at 31 evenly spaced spatial frequencies f between 0 and 0:4 mm −1 and three phases each separated by 2π=3.
Several experimental factors independent of the dynamics of the scattering particles can also effect the speckle contrast. These include speckle and pixel size, sampling window size, and shot noise [30]. In addition, finite source coherence length and polarization reduce the speckle contrast, which adds to the difficulty of obtaining quantitative measurements [31].
The total effect appears as a multiplicative correction factor for the contrast β. We adjust for this factor by first maximizing β to be as close as possible to unity. This was done by using a long coherence length (∼1 m) light source, and setting the speckle to pixel size at two pixels per speckle in order to meet Nyquist criteria. β was then found empirically by calculating the contrast from a rough metal object and found to be β ¼ 0:8. An additional reduction of 1= p 2 was also incorporated because of polarization loss [16].

B. Data Acquisition
The projected sinusoidal pattern through the LCOS can be written as where a DC and a AC are the weight of the planar offset and amplitude, respectively, of the projected wave. Next, we take advantage of the demodulation scheme proposed in [9]. This technique allows one to isolate the amplitude of a single frequency component by combining images at three equally separated phases φ 1;2;3 ¼ ð0; 2π=3; 4π=3Þ such that Next, we write the contrast remitted from the sample as where σ i and I i refer to the standard deviation and mean filtered image, respectively. K AC can then be extracted using the demodulation technique [9,32]: A summary of this data flow process can be seen in Fig. 2.

C. Liquid Phantoms
We measured three experimental systems. In the first, to help illustrate the necessity of incorporating optical properties into the speckle contrast model, three wells were filled with liquid phantoms using Intralipid fatty emulsions (Fresenius Kabi) as the scattering medium. Different concentrations of Intralipid and absorbing dye were added to achieve unique sets of scattering and absorption coefficients [33,34]: μ 0 s ¼ 1:00 mm −1 , μ a ¼ 0:00033 mm −1 ; μ 0 s ¼ 1:70 mm −1 , μ a ¼ 0:00033 mm −1 ; and μ 0 s ¼ 1:00 mm −1 , μ a ¼ 0:01 mm −1 . In this case, the Brownian diffusion of the scattering particles is constant, but the speckle contrast will be different in each well due to the change in optical properties.
To further validate our technique, liquid phantoms were created using 800 and 1026 nm polystyrene microspheres (Spherotech Inc.). The 800 nm spheres were diluted to 0.3% weight/volume (w/v) to achieve a reduced scattering coefficient μ 0 s ¼ 1:00 mm −1 at 633 nm wavelength, and Nigrosin dye was added to achieve an absorption coefficient μ a ¼ 0:005 mm −1 . Measurements of this phantom were done at exposures of T ¼ 10 ms. The 1026 nm spheres were diluted to a concentration of 1:1% w=v with corresponding μ 0 s ¼ 2:50 mm −1 , no absorber was added besides water, μ a ¼ 0:00033 mm −1 , and here T ¼ 1 ms. Microspheres are assumed to undergo diffusive dynamics, and the theoretical values of D b can be calculated using the Stokes-Einstein formula for a diffusing particle in a solution: where T is the temperature, η is the viscosity of the solvent, and r is the particle radius. The solution of microbeads is assumed to follow three-dimensional Brownian motion random walk statistics, where hΔr 2 ðτÞ ¼ 6D b τi, and D b is fit as the free parameter. Finally, liquid phantoms were also created using Intralipid and glycerin solutions, which allows one to adjust and control the viscosity of the solution. Intralipid was mixed with a water-glycerin solution at 1% concentration to achieve μ 0 s ¼ 1:00 mm −1 , and absorbing Nigrosin was added to achieve μ a ¼ 0:005 mm −1 . The glycerin to water ratios were then adjusted from 0%-43%.

RESULTS AND DISCUSSION
A. Monte Carlo Simulations Using our MC algorithm, PðY Þ was obtained for optical properties set at μ 0 s ¼ 1:00 mm −1 and μ a ¼ 0:005 mm −1 . The resulting distributions are plotted in Fig. 3 and were generated using 1 million photons, which took approximately 3 min on an Intel Core i7 processor. The distributions at short and long radial distance, ρ, illustrate the general shift of the mean momentum transfer quantity as one gets closer to the source. DWS begins to break down in this regime, as ρ→0 [21].
Inserting the above distribution into Eqs. (8) and (9) gives the autocorrelation as a function of radius and spatial frequency, respectively. This is compared to the correlation diffusion model from Eq. (7) for several moments in time τ, shown in Fig. 4. The autocorrelation has been normalized to unity at τ ¼ 0. As expected, the diffusion and MC models agree for short integration times and small k. For τ < 100 μs, the two models agree reasonably well for all simulated values of k. This is in agreement with the literature, which predicts the maximum valid time to be of the order of hundreds of microseconds [20]. However, LSI is usually performed with integration times of milliseconds, where the plots show significant divergence between the two models Thus MC modeling is necessary for quantitative analysis of LSI data. Fig. 2. Data flow for collecting AC speckle contrast amplitude. An image of actual data recorded at each step is provided for reference. The raw sinusoidal speckle image is collected at three equal phases with a CCD chip. Window mean and standard deviation filters of size 7 × 7 are then run over the images. The AC amplitude is extracted for each using the demodulation scheme. Finally, the contrast is determined by taking the ratio of the standard deviation to mean.

B. Optical Contrast
Intensity and speckle contrast maps for the three wells is shown in Figs. 5(a)-5(d) for T ¼ 10 ms. Increased absorption and scattering contrast can be seen in the top wells in the images. This is also shown in the plot of speckle contrast in Fig. 5(e). With a single scattering model that does not incorporate optical properties, the fitted D b varies between wells by up to 27% and are 1 order of magnitude greater than expected from the particle radius [34]. This is contrasted with the fit D b using our MC model, which has variation of less than 4% and agrees with literature values.

C. Polystyrene Beads
The resulting speckle contrast values are plotted as a function of spatial frequency alongside the MC and diffusion model fits in Fig. 6. The mean percentage error on the contrast fits for both particle radii are found to be approximately 4.55% and 0.86% for the diffusion and MC fit, respectively, showing better agreement with the MC model at this T value. This is due to the breakdown of the diffusion approximation assumption μ 0 s ≫ μ a;dyn in Eq. (7). In fact, for a 633 nm light source, T ¼ 1 ms, and D b ∼ 5:35 × 10 −7 mm 2 =s, μ a;dyn ∼ 0:6 mm −1 . The dynamic absorption here is actually of the order of the reduced scattering!    At room temperature T ¼ 298 K, η ¼ 9 × 10 −4 Pa⋅s, the predicted and fit diffusion coefficients for d 1 ¼ 800 nm and d 2 ¼ 1026 nm microspheres are shown in Table 1. The MC fits agree well with the theoretical values, and notice that the particle size ratios are preserved, in that This suggests the method has potential to provide accurate, quantitative assessments of particle motion in an imaging modality.

D. Glycerine Intralipid Titration
The viscosity η for glycerine and water solutions [35] are presented in Table 2. The expected fractional increase in diffusion coefficient is also shown, and is proportional to 1=η. These values are compared with the fit using MC modeling.
Here D 0 is the diffusion coefficient of Intralipid in pure water.
The fitted contrast is plotted as a function of spatial frequency in terms of D 0 in Fig. 7, and is found to fit well with a mean percentage error of 1.9%. Note that the predicted values with increasing viscosity are observed. The particle size of Intralipid has a large variance and depends on concentration and manufacturer [34,36], making it difficult to accurately predict the self-diffusion. It also contains various amounts of other solvents, which alter the inherent viscosity of the solution. Here D 0 is fit to ∼2 × 10 −6 mm 2 =s. This agrees with the general observation that the Intralipid particle radii are of the order of hundreds of nanometers [34], and smaller than the polystyrene beads from Subsection 4.C.
We believe the discrepancies in the fits, specifically around low frequencies, are likely due to the shape of the intensity profile. For most lasers, the beam profile is a Gaussian, which introduces additional frequency components into the image. In preliminary testing, we found that this can have a large effect, up to 15%, on the expected speckle contrast for planar illumination. We plan to conduct a more extensive study regarding the effect of the laser intensity profile on speckle contrast in the future, using concepts of cSFDI.

CONCLUSION
We have shown the capability for doing quantitative measurements of particle motion using the imaging technique LSI. MC simulations were carried out to obtain the photon autocorrelation as a function of spatial frequency, and compared to diffusion. Data were shown to agree with the values predicted using MC over the diffusion approximation. The Brownian diffusion coefficient was determined for polystyrene microspheres and shown to match theoretical values predicted by the Stokes-Einstein diffusion equation. Finally, by varying the viscosity of the scattering solution, the correct fractional increase in particle diffusion was observed.
In ongoing research, we aim to test this technique using a combined measurement system for reflectance and coherent SFDI that is currently being developed. We also wish to incorporate additional particle flow models and phantoms to simulate laminar dynamics. Static scattering components can cause measurements to become nonergodic, and we plan to employ multiple exposure techniques to separate static components from the image, if necessary [37,38]. We find that running a MC fit for each pixel individually is not feasible in terms of processing time, and we plan to use other developed methods, such as rapid-find lookup tables. Finally, we note that there are additional effects that spatial frequencies have on light propagation that may be advantageous. In particular, it has been shown that penetration depth is also a function of spatial frequency [9,23]. Thus, tomography may be performed using SFDI [11] to reconstruct three-dimensional maps of speckle contrast.