Determination of the effect of source intensity profile on speckle contrast using coherent spatial frequency domain imaging.

Laser Speckle Imaging (LSI) is fast, noninvasive technique to image particle dynamics in scattering media such as biological tissue. While LSI measurements are independent of the overall intensity of the laser source, we find that spatial variations in the laser source profile can impact measured flow rates. This occurs due to differences in average photon path length across the profile, and is of significant concern because all lasers have some degree of natural Gaussian profile in addition to artifacts potentially caused by projecting optics. Two in vivo measurement are performed to show that flow rates differ based on location with respect to the beam profile. A quantitative analysis is then done through a speckle contrast forward model generated within a coherent Spatial Frequency Domain Imaging (cSFDI) formalism. The model predicts remitted speckle contrast as a function of spatial frequency, optical properties, and scattering dynamics. Comparison with experimental speckle contrast images were done using liquid phantoms with known optical properties for three common beam shapes. cSFDI is found to accurately predict speckle contrast for all beam shapes to within 5% root mean square error. Suggestions for improving beam homogeneity are given, including a widening of the natural beam Gaussian, proper diffusing glass spreading, and flat top shaping using microlens arrays.

1. Introduction

Dynamic Light Scattering (DLS) is a method which utilizes the decorrelation of coherent light to estimate the motion of particles in scattering media. DLS is used across scientific disciplines to characterize samples in physical, chemical, and biological settings [1-3]. In biomedical applications, several DLS methods are used, including: Laser Doppler flowmetry [4-6], Diffusing Wave Spectroscopy (DWS) [7-10], and Laser Speckle Imaging (LSI) [11-15]. Common to all techniques is that decay of the autocorrelation function of coherent light is related to the degree of scattering particle motion.

It has been shown that in highly scattering media, the electric-field temporal autocorrelation function G1(τ) = 〈E(0)E*(τ)〉 obeys a transport equation [16]. Using this model, the autocorrelation function can be calculated and compared to experimental measurements in order to determine the dynamics of the light scattering structures. This approach has been used with simple, point like sources extensively for fiber based techniques such as DWS, in which the distance between the light source and detector is clearly defined. In contrast, for imaging based measurements such as LSI, source-detector separations are ill defined. However, by taking advantage of the spatial filtering characteristics that turbid media enact onto structured light, one can model correlation transport as a function of the spatial frequencies in the structured light projected on to the sample.

This concept has been applied to the transport of light radiance as the foundation for a reflectance imaging technique called Spatial Frequency Domain Imaging (SFDI), which has been well established and used to determine the absorption, scattering, and fluorescent properties in turbid media [17-23]. Recently, we applied experimental methods and modeling precedent from SFDI to the transport of correlation by projecting structured coherent light on dynamic turbid samples, and measuring the speckle contrast as a function of projected spatial frequency [24]. We showed that the speckle contrast follows expected behavior as a function of spatial frequency [25], and that correctly modeling the effect using correlation transport allows one to fit quantitatively for the dynamical properties of light scattering particles [24].

The dependence of speckle contrast on spatial frequency shows that despite the fact LSI may be independent of overall intensity of the laser source, spatial variations across the beam spot itself can cause the speckle contrast to change. In typical LSI experiments the beam has a naturally Gaussian intensity profile, is projected over a finite extent, and may have additional artifacts from traveling through a beam spreader such as a lens or glass diffuser. Thus, to increase the quality of speckle imaging experiments, it is important to study the nature and magnitude of the phenomenon more in-depth. This manuscript first demonstrates effects of the beam profile qualitatively using in vivo experiments, then performs quantitative predictions and analysis using coherent Spatial Frequency Domain Imaging (cSFDI) of controlled liquid phantoms.

2. Theory

The speckle contrast parameter is defined as

where K is the speckle contrast, σ is the intensity standard deviation, and 〈I〉 is the mean intensity. When imaging the speckle pattern with some finite integration time T, K is related to the electric field autocorrelation function G1(τ) through the expression [14,26]

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. If the scattering objects are static, the field does not change and K reduces to β. Assuming the scatterers are in motion, a correlation diffusion equation can be used to find G1 analytically as a function of optical properties and the mean square displacement 〈Δr2(τ)〉 [10]. However, we and others have found 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 [24,27]. Therefore, we utilize Monte Carlo (MC) simulation techniques to determine the autocorrelation function. As proposed in [9,10,28], 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 [29]:

Here P(Y) is the distribution of dimensionless momentum transfer, Y = 1 – cos(θ), and k0 is the wave number inside the medium.

To take advantage of the cSFDI technique, the autocorrelation is then calculated as a function of spatial frequency by applying a numerical Hankel Transform to G1 [24]

Finally, the speckle contrast modulation transfer function K(f), can be found by inserting G1 into Eq. (2). K(f) is a spatial filter that the sample enacts onto the speckle contrast image when the coherent source has spatial structure.

3. Methods

Unless otherwise noted all experiments were performed using a 30 mw 633 nm HeNe Laser (JDS Uniphase, Inc) coherent light source, and images were captured with a Retiga EXI cooled 12-bit CCD camera (Qimaging, Inc) at an integration time of T = 15 ms. The optical projection 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. Speckle contrast was calculated using a sliding window mean and standard deviation filter, calculated over small sampling neighborhoods of 7x7 pixels. Border areas within 7 pixels of the edge were discarded. Speckle flow index was defined as the reciprocal contrast FI = 1/K. All measurements were carefully designed to meet or exceed Nyquist criteria to minimize sampling errors [30]. The value of β was found empirically using a static diffusive target and calculated to be β = 0.8. Note that finding β in this manner may also incorporate shot, readout, and other types of noise which can inflate the value. However, with adequate signal and cooling this effect should be small, and any error in calculating β will appear as a multiplicative factor in images of speckle contrast and not alter the conclusions of the paper at large.

Two in vivo speckle contrast measurements were performed. In one experiment, images were acquired through the skin on a human radial artery, near the superficial palmar branch, with a Gaussian illumination profile of FWHM 10 mm. To track location, opaque fiducial markers approximately 2 mm in diameter were placed near the artery using adhesive. The branch point of the artery was first centered under the peak of the Gaussian, then moved diagonally outward towards the tail edge of the beam. A second in vivo measurement was taken on a mouse cortex, where Gaussian and flat-top laser beam profiles were used to image a 17-month old C57/BI6 male mouse under 1% isoflurane anesthesia. A bilaterial craniectomy was performed to expose the brain.

Experiments were also performed with liquid scattering phantoms of known optical and dynamical properties composed of a microbead solution. Phantoms were mixed to a concentration of 0.48% w/v with 800 nm Polystyrene beads (Spherotech Inc.) and distilled water. The optical properties, calculated using standard Mie Theory and known water spectra Motion diffusion characteristics, where the mean square displacement 〈Δr(τ)〉 = 6D and D can be calculated using the Stokes-Einstein formula; in this case D6.1x10− 7 mm2/s. Speckle contrast data was taken using the following beam shapes: standard Gaussian profile, a flat-top profile shaped using a fly’s eye lens array pulled from a digital projector unit (Epson, Inc), and the output through a frosted glass diffusing plate (Thorlabs, Inc). Note that for coherent light the intensity pattern formed through a diffuser is speckle in nature and depends on the radius of the beam waist and the distance between the diffuser and object [26]. The beam width used here was typical for a HeNe output, with 0.5 mm FWHM, and the diffuser placed at a distance of 10 cm from the sample. Control measurements were taken with each beam shape on a Spectralon reflectance standard. This was done to confirm that variation in speckle contrast is not caused by systematic errors such as low light intensity, noise, or vignetting. Spectralon is composed of static highly scattering material where no change in speckle contrast is expected with respect to beam shape.

To numerically simulate cSFDI data, a simple time resolved Monte Carlo photon propagation algorithm was obtained from the Oregon Medical Laser Center Website [31], and a new feature was added to track the weighted dimensionless momentum transfer. Ten million photons were propagated through a semi-infinite medium with optical parameters matching our liquid phantom above, and P(Y) found at 200 discrete radial bins linearly spaced between 0 mm and 40 mm. G1(f) was then calculated using Eq. (3) and Eq. (4), and finally the transfer function of the contrast image K(f) calculated with Eq. (2). A plot of K(f) is shown in Fig. 1 for a typical reduced scattering coefficient in tissue μ1.7 mm1. Note that the scattering medium acts as a high pass filter on the speckle contrast.

Fig. 1

Plot of speckle contrast as a function of spatial frequency using Monte Carlo simulations. This is the “speckle contrast modulation transfer function” K(f), which is applied as a spectral filter on beam profile images to predict remitted speckle contrast.

The Fourier spectral coefficients of the reflectance image for each beam shape were determined by implementing a Discrete Fourier Transform algorithm. For finite exposures, the recorded speckle is a stochastic process and thus the mean-filtered image 〈I〉 has nonzero variance which depends inversely on the number of sampled pixels, determined by window size, according to the central limit theorem. This variance manifests itself in the high spatial frequencies when using DFT. To decouple the intensity profile’s spatial frequencies from this natural statistical variance, over 100 independent images of 〈I〉 were averaged together, effectively increasing the sample number by two orders of magnitude. The diffuse intensity image was then represented with a Fourier series,

The same spectral representation was then applied to the standard deviation filtered images, except the speckle contrast MTF served to reduce each component. Rearranging Eq. (1),

The predicted speckle contrast image is then just the ratio: K〈I〉.

4. Results

Figure 2(a) shows a human radial artery illuminated under different sections of a Gaussian beam. Raw images and generated flow maps are shown in Figs. 2(b) and 2(c).

Fig. 2

a) Human radial artery illuminated with a Gaussian beam. b) Arterial branch placed directly beneath the peak (left), then moved diagonally and placed underneath the Gaussian tail (right). Dark ellipses represent fiducial markers. c) Visual differences in flow index are seen between the two locations. These depend on region of interest, but can reach 13% in the high flow arterial sections. d) A profile plot taken horizontally across the image starting at the vertical center further illustrates differences.

One observes a bulk increase in calculated flow index of approximately 15% as the vessel branch is moved radially outward, away from the peak of the beam profile, as the photons exiting the tissue here have travelled longer average path lengths. This occurs because the measurement point is made further from the illumination peak at the center of the beam.

This effect was studied more closely using homogenous liquid phantoms, where DFT’s could be used to generate a forward contrast image using Monte Carlo simulations of cSFDI. Results are showed in Fig. 3 for three common beam profiles, a standard Gaussian shape, output through a diffusing glass plate, and a flat homogenous profile. The reflectance averaged over 500 exposures is shown in Fig. 3(a). Figure 3(b) show the experimental speckle contrast calculated with the usual sliding window standard deviation and mean filters. These are to be compared with the predicted values using cSFDI (Eq. (5) and Eq. (6)) shown in Fig. 3(c).

Fig. 3

(a) Intensity image for three distinct beam shapes, a Gaussian shape, output through a frosted glass diffuser, and flat top. (b) Experimental speckle contrast images computed using 7x7 neighborhoods. (c) Predicted speckle contrast using the forward cSFDI model based on the spatial frequencies of the beam shape. Note that artifacts within the beam shape are amplified due to the high pass nature of the speckle MTF. (d) (e) Effects from systematic variables such as vignetting, sliding window filter response, or low counts are shown to be negligible by direct comparison with a static object, where the speckle contrast is seen to be constant and approximately equal to β, as expected.

It is clear from the images that cSFDI is able to reasonably predict the behavior of the speckle contrast. The root mean square error (RMSE) was calculated at each pixel and normalized by the theoretical value predicted by cSFDI to obtain a percent error. The average percent error over all pixels was less than 5% for all three shapes.

To act as a control, speckle images were also taken of static scattering objects using a Spectralon reflectance standard. The raw images and calculated contrast are shown in Figs. 3(d) and 3(e). As expected, the calculated contrast is effectively uniform across the image space regardless of beam shape and approximately equal to β.

Figure 4(a) shows the raw speckle image of a mouse brain for two intensity profiles, a Gaussian and flat shape. Speckle flow maps were generated and are shown in Fig. 4(b). Although the large variations within the image between vessels and tissue conceal the effect, inspection of a single section of the superior sagittal sinus vessel in Fig. 4(c) shows a bulk change of approximately 10% in calculated flow index between the two images. Finally, the percent difference image between the two flow maps is calculated across the cortical surface and shown in Fig. 4(d).

Fig. 4

Images of mouse brain with a Gaussian beam shape (left) and flattened beam shape (right). (a) Raw images are shown where the shape of the beam is evident. (b) Flow index maps and (c) a specific ROI of super saggital sinus vessel illustrate a difference in perceived venous flow that depends on the beam shape. (d) The percentage difference between flow maps show significant variation across the field.

Here, the differences in flow caused by the beam shape are clear, showing changes in some parts by as much as −25% near the peak of the Gaussian, due to shorter average path lengths, and + 25% looking further out along the tail. This agrees with predictions made by cSFDI and the in vivo human artery experiment presented above.

5. Discussion and conclusion

Our results show that the source intensity profile can have a significant effect on the speckle contrast and thus on the resulting speckle flow index. This phenomenon occurs because of variations in the average path length of photons across the beam profile. When measuring near a local intensity maximum, photons generally need to travel shorter path lengths to reach the detector, and thus transfer less overall momentum, as compared to lower intensity regions. The effect causes changes in the measured flow index across the field of view, and if not accounted for or minimized, would confound any attempt to determine absolute rates of blood flow.

Under controlled conditions, the close agreement between the forward model predicted by cSFDI and experimental results suggest this technique is an accurate method for predicting the effect of light structure on speckle contrast. The relative ease with which the speckle contrast filter K(f) can be applied to any complicated structure using DFT’s, such as the “diffuser shape” in Fig. 3, illustrate the utility of the method in an imaging geometry. In principle cSFDI could be used to generate a multiplicative “correction” factor to account for beam shape. However, this creates a number of extra complications when imaging a heterogeneous sample. First, diffuse beam shape can no longer be analyzed using DFTs, as heterogeneities would affect the Fourier spectrum. In addition, the optical properties must be known at each location in the image and the speckle MTF and forward model must generated for each set. If the beam shape and bulk optical properties were assumed a priori, a simple correction factor could be generated, but this approach would contain errors that scale with the magnitude of the heterogeneities. Perhaps the ideal solution is to remove structure from the beam shape before measurements in the interest of quality data collection. Applying this philosophy to an in vivo experiment in a mouse brain shows major differences when taking steps to appropriately homogenize the beam.

There are several methods to help reduce variation in the beam profile. The simplest is to expand the beam to a larger radius, effectively increasing the width of the natural Gaussian profile. The disadvantage to this method is a significant decrease in the intensity per unit area, and that the addition of optical elements can create more artifacts and aberrations. Another possibility is expanding the beam radius before entering a diffuser plate and decreasing the distance to the object. It is known that coherent light through a diffuser will project an intensity speckle pattern on an object some distance away, where the speckle size is governed by the incident spot radius and the projection distance. If this secondary intensity speckle pattern is fine enough such that the Nyquist criterion is not met, the path length distribution will be uniform across the pixel array due to under sampling. For a circular spot of diameter D incident on a diffuser and projected a distance z away, the power spectral density is defined [26] to be a bandwidth limited function with maximum frequency f:

Thus, to avoid sampling the intensity speckle pattern with the camera array, the pixel width x must be larger than the Nyquist limit:

Note this under sampling serves to blur the projected secondary speckle in conjunction with standard light diffusion through tissue, creating a very smooth profile. However, while the Gaussian profile of the beam is broadened when passing through a diffuser, it is not eliminated.

A third option is to use a flat top beam shaper with a homogenizing element. Commercial ones are typically manufactured using a fly’s eye micro-lens array, and these optical systems are also used in commercial projectors to form a homogenous rectangular beam. Differences in path length between elements of the lens array can cause a secondary interference speckle-like pattern [32], but our experience has been that these have been under sampled and blurred, as in Eq. (8), to create the homogenous flat shape seen in Fig. 3.

In conclusion, the shape of the source beam causes a significant effect on the speckle contrast and speckle flow index measured in LSI. We demonstrated this effect using two types of in vivo measurements: small animal brain, and human radial artery. Using controlled microsphere phantoms, the effect was quantified and predicted using cSFDI framework. This was done both for a Gaussian beam and for the more complicated intensity profile from a diffusing glass plate. We showed that spatial frequency domain analysis provides a means to predict the effect of any beam shape on the speckle contrast by calculating a speckle contrast MTF. We also investigated and reported methods for increasing the beam homogeneity in LSI experiments.