Laser Speckle Modeling and Simulation for Biophysical Dynamics: Influence of Sample Statistics

. Laser speckle is a statistical phenomenon and should be treated as such when establishing experimental parameters for using laser speckle techniques to infer information about the dynamics of a coherently illuminated medium, such as biological tissue. Herein, we demonstrate that when sampling (imaging) speckle patterns, it is critical that the sampled speckle patterns are unbiased estimators of the underlying random speckle field. If this condition is not met, the quality of the results will be compromised. Specifically, this study examines the effects of first and second order spatial statistics of speckle intensity images on laser speckle contrast imaging results. Finally, it is recommended that when using speckle techniques such as laser speckle contrast imaging, investigators should examine the first and second order spatial statistics of a speckle image prior to collecting actual data. If this examination reveals that the imaged speckle intensity image is not an unbiased estimator of the underlying random speckle field, then adjustments should be made to ensure that the images taken are unbiased estimators of the true speckle field. © 2017 Journal of Biomedical Photonics & Engineering.


Introduction
Laser speckle techniques are well established for investigating biophysical dynamics, including the motion of blood cells [1][2][3] and tissue mechanics [4][5]. Laser speckle techniques possess numerous advantages as tools for investigating dynamic biological systems, including but not limited to good spatial resolution, high sensitivity to motions, relative simplicity, and cost effectiveness. These same features, however, present a weakness of the techniques. That is, it is relatively simple to get data. The simple act of illuminating a scattering surface or volume with coherent radiation (i.e., laser light) will result in the appearance of a speckle pattern in the observation plane. If the scattering medium is dynamic, then the observed speckle pattern will also appear to be temporally dynamic and estimates of this motion are easily made [6].
However, laser speckle is a statistical phenomenon and this fact must be accounted for in both the experimental parameters that are established to spatially and temporally sample the scattered speckle field and also in the approaches used to assess the speckle dynamics. In other words, the sample statistics of the observed intensity pattern must be unbiased estimators of the true scattered field. The result of not meeting this minimum sampling requirement can result in significant errors in the analysis of speckle dynamics. This paper will focus on the spatial sampling of the speckle intensity patterns and will not address temporal statistics.
The first and second order statistics of laser speckle patterns are well established for both 'fully developed' (polarized) and 'partially developed' (un-polarized) speckle patterns [7] and it is not necessary to present a full description of these statistics herein. Yet, a brief overview is necessary as the remainder of the present paper explores the effects of inadequate spatial sampling of laser speckle patterns such that the sample statistics are not unbiased estimators of the true underlying statistics.
We begin by noting that laser speckle may be viewed as the result of the coherent summation of randomly phased Huygen's wavelets. Indeed, speckle as a general phenomenon arising from the coherent illumination of a scattering surface or volume can be viewed in this fashion regardless of the form of the radiation (e.g., coherent x-ray, ultrasound, etc.). Thus, for the remainder of this paper we will generalize and simply refer to the speckle phenomenon as speckle and drop the 'laser' modifier. It is left up to the reader to make the appropriate substitutions when necessary. These randomly phased wavelets can be viewed as arising from a real source in the case of imaged (or subjective) speckle or a virtual source in the case of non-imaged (or objective) speckle. Observing the speckle pattern intensity with a bare CCD or CMOS chip without an intervening system of imaging lenses would be an example of the latter case. Note that this configuration is equivalent to imaging at infinity. Regardless of the mode of observation, however, the physics of the speckle formation is identical and both imaged and non-imaged speckle display identical first and second order statistics. For this manuscript, first order statistics refer to the statistics at a single point in time and space, while second order statistics describe the relationships between neighboring points. First order statistics describe, for example, the intensity at a given pixel, while second order statistics describe the spatial structure of a speckle pattern [8].

First Order Statistics
Consider the situation displayed in Fig. 1. The desire is to describe the fields in some observation plane L distant from the illuminated surface. To do this, we consider the surface fields to be equivalent to a finite number of secondary Huygens' wavelets. The 'finite' restriction will be subsequently relaxed. If we ignore the n R / 1 dependence of the individual wavelets and make the valid assumption that there is no systematic J of Biomedical Photonics & Eng 3 (4) 31 Dec 2017 © JBPE 040302-3 variation in phases (i.e., the wavelets are truly randomly phased) then the fields at the surface can be expressed by where N is the number of summed wavelets. Fig. 1 Geometry of speckle observation. The observation plane is located at L and s is the extent of the illuminated portion of the surface.
Since we are treating the observation fields as the sum of randomly phased phasors, as N → ∞ the sum becomes asymptotically circularly Gaussian. The joint probability density function (PDF) for the real and imaginary parts of a circular Gaussian complex variable is readily shown to be [7,8] where σ 2 = σ r 2 = σ i 2 is the standard deviation of the underlying random process. By integrating over the real and imaginary parts such that the integral of the joint PDF equals unity and switching to polar coordinates where The term in the brackets of Eq. (3) is the PDF of the amplitude and defines the Rayleigh distribution when a ≤ 0 .
Since intensity, not amplitude, is what is observed in Fig. 1, by following the preceding approach for intensity, I = a 2 , the integral becomes The term in the brackets is the negative exponential PDF of intensity (Fig. 2). This PDF can be expressed by [9] p I i

Laser Speckle Contrast
The contrast, c, of a laser speckle pattern is a first order property that is frequently investigated in order to infer properties related to the dynamics of a biophysical system. For example, laser speckle contrast imaging (LSCI) is a commonly employed modality used for the relative and qualitative imaging of blood flux [1,2,10]. LSCI and related modalities are relatively simple to implement, have wide fields of view, and relatively good spatial and temporal resolution. These features have allowed these techniques to be used as powerful tools in the measurement and monitoring of biophysical systems in near real-time.
In addition to the statistical considerations discussed below, LSCI results depend upon several other experimental considerations as discussed in the review  by Briers et al. [2]. Recently, Khaksari and Kirkpatrick [10,11] have demonstrated that ultimately, LSCI is sensitive to advective flux, scattering and absorption. Indeed, these three variables are challenging to disassociate from one another in terms of their influence on LSCI.
The contrast, c, of a speckle pattern is typically defined as the quotient of the standard deviation σ I and the mean µ I of the measured intensity: For speckle patterns exhibiting and exponentially distributed intensity PDF (see Fig. 2), this value is exactly unity. The concept behind LSCI is to relate the motion of scatterers (blood cells) in the object space to the motion (and thus time integration) of speckles in the observation, or imaging space. Thus, the imaged speckle intensity pattern must be an unbiased estimator of the true dynamic speckle field arising from the moving scatterers. Otherwise, establishing this relationship between the moving scatterers and the observed speckle intensity is fundamentally impossible.
Thus, while contrast is defined in terms of expected values (see Eq. 6), in practice it is calculated over a local spatial (or temporal, or both) neighborhood in terms of the sampled statistics for the intensity: where M is the mean intensity of the sampled speckle and S 2 is the standard deviation of the sampled intensity, I. If it is assumed that the sample statistics (M and S 2 ) are unbiased estimators of the true field, then as the number of samples theoretically goes to infinity, then the estimate of the contrast approaches the theoretical value of unity. In practice, the number of samples does not need to approach infinity, it is sufficient at this point to note that the number of samples just needs to be 'large'. Thus, we can write: This is demonstrated in Fig. 3. The details of how the speckle pattern in Fig. 3 was generated are given below in the Methods section, as is the calculation of minimum speckle size.
As demonstrated by Duncan et al. [12], and is obvious from the contrast image in Fig. 3, the local values of contrast exhibit a distribution. In Fig. 3, the local contrast (calculated over a 7 × 7 sliding neighborhood window) values ranged from a minimum of approximately 0.2 to a maximum of approximately 1.7. An examination of the PDF of the contrast values reveals that the PDF of the local contrast is closely described by a log-normal distribution of the form where c m is the median value and σ g is a dimensionless width parameter. The shape of this PDF is shown in Fig. 4   It should be noted here, that typically when discussing sample statistics, the implicit assumption (although not required) is that the samples are statistically independent. Recall that in the example above, the speckles size equaled 4 pixels. This is 2x the spatial Nyquist criteria [13]. Thus in the example above, the speckle field was sampled above the spatial Nyquist limit and the smallest speckle extended over several pixels. Also note, that as is discussed below, the correlation length of the speckle pattern is on the order of the smallest speckle. Thus, the samples of intensity are not statistically independent.
This lack of statistical independence suggests that the median c m and width σ g parameters of Eq. 9 will depend upon the size of the neighborhood over which the contrast is computed and on the size of the speckles relative to size of the pixels in the observing camera. Absolute size of the speckles is not a consideration in and of itself, per se. It is the relative size of the speckles to the pixels that is of interest. These considerations are examined below.

Second Order Statistics
The spatial structure of a speckle pattern is described by the size of the observed speckles. Speckle size can be approximated as the square root of the coherence area of the speckle intensity pattern. The coherence area, in turn, is a function of the covariance of the intensity distribution [8]. The second order statistics of a speckle pattern describe, in essence, the coarseness [7] of the spatial structure. The distribution of speckle sizes in a speckle pattern has readily been found by either examining the full width at half-maximum (FWHM) of the autocorrelation function of the speckle intensity pattern or, alternatively, the spatial frequency distribution of the speckle pattern can be quantified by calculating the power spectral density (PSD) of the intensity pattern [7,14]. Either approach yields the same information as the autocorrelation function and the PWD constitute a Fourier transform pair. An alternative approach of using spatial Poincaré plots with variable lag lengths was recently presented by Majumdar and Kirkpatrick [15]. This approach reduces some of the uncertainties associated with the approaches discussed above, particularly as the speckle size gets large and also provides a more physical connection with speckle size and the correlation length of the intensity pattern. The effects of spatially under sampling the speckle filed is discussed in subsequent sections of this paper.
The remainder of this paper is devoted to a discussion of the effects of spatial sampling of laser speckle fields on laser speckle contrast images. In particular, the remainder of the manuscript focuses on the situations where the sample speckle image is not a statistically unbiased estimate of the true underlying scattered field. Both simulated speckle fields and actual, imaged speckle intensity patterns are used as examples.

Laser Speckle Numerical Simulation
Simulated band-limited speckle patterns were generated in the Matlab environment by filling a circular masks of radius R r embedded in a larger zero-padded space with random complex numbers over the interval of 0,2π ⎡ ⎣ ⎤ ⎦ .
The spaces were Fourier transformed and multiplied by their complex conjugates, generating polarized speckle patterns exhibiting negative intensity PDFs (Fig. 2). The size of the speckles was altered by changing the relative size of the circular mask to the surrounding zero-padded space. When the speckle size equaled 2 pixels, the spatial Nyquist criteria was met [13]. This process for generating synthetic speckle patterns has been used extensively in the past [2, 6, 12, 13-15]. Using these speckle patterns, the effects of the neighborhood size over which speckle contrast was calculated was investigated as a function of speckle size and the results compared to those obtained when the sampled intensities were statistically independent from each other. That is, when the speckle size ≤ 1.0 pixels (see above). Neighborhood sizes ranging from 3 × 3 pixels up to 13 ×13 pixels were investigated. The largest speckle size investigated was 3 × the Nyquist criteria size (i.e., 6 pixles) and the smallest (undersampled) speckles were 1 pixel in size.
These simulations were used explicitly to examine what happens when the second order statistic (speckle size) of the sampled speckle intensity pattern is not an unbiased estimator of the underlying speckle field.

Laser Speckle Contrast Imaging Simulation Using a Spatial Light Modulator
The effects of incorrect first order statistics on LSCI was investigated by simulating LSCI experiments with a nematic liquid crystal spatial light modulator (SLM) operated in the phase-only mode [16]. The focus of these experiments was to adjust the total intensity of the speckle patterns created via the SLM in the image plane of a CMOS camera and assess the resulting first order statistics of the speckle patterns. As the intensity PDFs of the speckle patterns changed as a result of altering incident intensity on the SLM, the effects on LSCI were investigated. Two separate analyses were performed to investigate the influence of first order statistics on LSCI. First, the global contrast of the speckle patterns generated with varying incident intensities was calculated in order to investigate whether the first order statistics had an influence on global contrast. Second, an analysis was performed on speckle patterns generated with the SLM that had a central region with a In the SLM arrangement, the SLM (Boulder Nonlinear Systems XY Series, model P512) is illuminated by a 2.5 mW stabilized HeNe laser emitting at 633 nm. The resulting speckle patterns were imaged via a 4f lens system with an 8-bit CMOS camera (Thorlabs, DCC1545M High Resolution USB 2.0 CMOS Camera). The half-wave plate ( λ / 2 plate) in the incident arm combined with the analyzer in the imaging arm allowed the SLM to operate in the phase only mode and also to produce polarized speckle patterns.
The magnification of the system was calculated and measured to be 0.5.
To generate speckle patterns from the SLM, a sequence of 50 phase screens was loaded to the SLM. The phase screens were generated in the Matlab environment by filing a 512 × 512 array with randomly generated complex numbers of unit amplitude with phases uniformly distributed over the interval 0,2π ⎡ ⎣ ⎤ ⎦ .
A look-up table ensured that the phase modulation resulted in a linear relationship between the phase distribution, 0,2π ⎡ ⎣ ⎤ ⎦ , and the resulting gray values 0,255 ⎡ ⎣ ⎤ ⎦ . The scattered light remained both temporally and spatially coherent, allowing for random interference to occur when the scattered light impinged upon the CMOS chip. Speckle patterns were thus produced. Using the copula approach described by Duncan and Kirkpatrick [14], the phase screens, and thus the resulting speckle patterns, de-correlated in a userdefined fashion. In the present experiments the speckles decorrelated following a Gaussian decorrelation function. Regions that were made to simulate a blood vessel decorrelated approximately 16 times faster than the surrounding, more slowly moving speckles. This value was chosen through a purely empirical process so that the resulting laser speckle contrast images visually resembled actual speckle contrast images acquired in J of Biomedical Photonics & Eng 3 (4) 31 Dec 2017 © JBPE 040302-7 Fig. 6 Dependence on local neighborhood size for the first two moments of the local contrast. Median (left) and width parameter (right) for a sampled speckle pattern exhibiting an exponential PDF. Speckle size was 4 pixels/speckle. Fig. 7 Parameters of contrast distribution as a function of local neighborhood and speckle sizes. The results for the statistically independent situation is presented as well (dotted line). the laboratory. The minimum speckle size in the imaged speckle patterns was approximately 3.25 pixels per speckle. Thus, the speckle size exceeded the minimum as determined by the (spatial) Nyquist criteria.
To create time integrated speckle patterns in order to create laser speckle contrast images, sequential frames were summed on an intensity basis to simulate the speckle blurring observed in an actual LSCI experiment when dynamic speckle imaged with a finite exposure time. The blurred speckle resulted in a reduced contrast as is seen in LSCI.
As the goal in these experiments was to alter the first order statistics of the sampled (imaged) speckle patterns, the incident intensity was varied over 4 orders of magnitude. This resulted in some speckle images with an excess of dark pixels (due to under-exposures) and at the other end of the spectrum, some speckle images with an excess of saturated pixels. Both of these conditions alter the first order statistics of the speckle images.

Effects of Neighborhood Size on LSCI -Second Order Statistics
To examine the effects of second order statistics on the PDF of a speckle contrast image, we first present the results of results of local neighborhood size over which speckle contrast is calculated [12]. As discussed above, the initial speckle size in this analysis was 4 pixels per speckle, or 2 × the Nyquist lower limit. As can be seen in Fig. 6 To examine the effects of speckle size (which is a second order statistic) on the above results, the simulation was repeated for varying speckle sizes. The results are presented in Fig. 7 for speckle sizes at the Nyquist limit, at 2 × and 3 × the Nyquist limit.
Additionally, the simulation was repeated when the sampled speckles were statistically independent of one another. In other words when the speckles were undersampled and had sizes ≤ 1.0 pixel.
From a practical standpoint, these results indicate that in order to maintain statistical validity, the speckle pattern must be sampled spatially in a manner that meets the Nyquist criteria. That is, the minimum speckle size must be at least 2 pixels per speckle.

Effects of Intensity PDF on LSCI -First
Order Statistics. Global Contrast As discussed above, global contrast on a properly sampled speckle image should equal unity. In this section, the results of sampling the speckle field so that the intensity distribution is not an unbiased estimator of the true scattered field is examined. As discussed above, this was accomplished by either under or over exposing the speckle images. The top set of images in Figs. 8 were taken at 100% power from the laser and the bottom set of images were taken with a relative intensity of 0.1× that maximum intensity (i.e., 4 orders of magnitude lower).
When the global contrast of speckle images with different exposures was calculated, it was found that those speckle images that had displayed a negative exponential intensity PDF had the highest contrast, while those that were either over exposed or under exposed, and thus did not display a negative exponential PDF, had a lower global contrast. These results are shown in Fig. 9.
Recalling that contrast is defined as the ration of the mean of the intensity to the standard deviation of the intensity (Eq. 6), then the results displayed in Fig. 9 can be explained by looking at each of these statistical descriptors individually. Figure 10 plots the means and standard deviations of intensity for the speckle images used to generate the data of Fig. 9. At the upper and lower ends of the plot, it is seen that the standard deviations decrease relative to the means as a result of pixel saturation (for over exposed images) or under exposure at the opposite end.   Likewise, for the simulated laser speckle contrast images, incident intensity that moved the intensity PDF away from a negative exponential PDF, resulted in a smaller difference between the 'fast' decorrelating area that mimicked a blood vessel and the slower decorrelating areas that mimicked the surrounding tissues. This is shown in Figs. 11a and 11b. Fig. 11a Simulated LSCI image. The central 'flow' region can be seen and the had a decorrelation time approximately 16 times faster than the surrounding 'static' region. Fig. 11b The ability to discriminate between the flow and static regions in Fig. 11a decreased when the incident intensity was either too low or too high. Red stars display the contrast in the flow region, while the black dots display the contrast of the 'static' region. In the central region of the curves, the intensity PDFs of the speckle images followed a negative exponential behavior. At the ends the intensity PDFs were not negative exponentials. Reproduced with permission from Ref. [16].
When the ratio of the contrast in the static region to the contrast in the flow region was plotted as a function of incident intensity, the trend seen in the above simulations was repeated. That is, when the intensity PDF of the speckle images was a negative exponential, the ratio in contrast between the two regions was a maximum. This is shown in Fig. 12. The ratio of contrast in the two regions was defined in decibels as Fig. 12 Ratio between static and flow regions in the simulations as a function of relative incident intensity.
As the ends of the curve, the ratio decreases. In the middle of the curve, the intensity PDFs of the speckle patterns were negative exponentials.

Discussion
The results of these simulations and experiments reinforce the fact that laser speckle is a statistical phenomenon and must be treated as such. The results further demonstrate that the sample speckle intensity must be an unbiased estimate of the true, underlying random process that speckle is. The results presented herein demonstrate that contrast is maximized the speckle field is sampled so that the sampled speckle intensity is described by a negative exponential PDF and so that the minimum speckle size is at least 2 pixels per speckle in order to satisfy the Nyquist criteria. Frequently, the use of the so-called 'coherence' factor, β , has been invoked in LSCI for the purpose of correcting for insufficient spatial sampling of the field [1,17,18]. This factor originates from the Seigert relation that relates the intensity autocorrelation function to the field autocorrelation function in quasielastic light scattering (QLS that it is statistically more appropriate to sample the underlying speckle field adequately in the first place, than to attempt to correct for this with a correction factor. The results of this study lead to a recommendation that prior to taking speckle measurements, such as LSCI measurements, that a single speckle image (single frame) be acquired and the first and second order statistics analyzed to ensure that the sampled speckle pattern is an unbiased estimator of the true, underlying speckle field. If the acquired speckles are too small, or if the intensity distribution of the speckle image is anything but a negative exponential (for a polarized speckle pattern), then adjustments should be made to ensure that these conditions are met prior to acquiring data.

Disclosures
The authors declare that there are no conflicts of interest related to this article. Portions of this paper were presented as part of the 2017 Saratov Fall Meeting, Saratov Russia, September 26-29, 2017.