Simulation of speckle patterns with predefined correlation distributions

: We put forward a method to easily generate a single or a sequence of fully developed speckle patterns with pre-defined correlation distribution by utilizing the principle of coherent imaging. The few-to-one mapping between the input correlation matrix and the correlation distribution between simulated speckle patterns is realized and there is a simple square relationship between the values of these two correlation coefficient sets. This method is demonstrated both theoretically and experimentally. The square relationship enables easy conversion from any desired correlation distribution. Since the input correlation distribution can be defined by a digital matrix or a gray-scale image acquired experimentally, this method provides a convenient way to simulate real speckle-related experiments and to evaluate data processing techniques.


Introduction
Speckle phenomena exist in all coherent imaging and non-imaging systems. Besides being experimentally investigated, simulation can be a powerful complement for understanding the properties of speckles, for evaluating data processing methods and for investigating the possibility of further applications of laser speckle phenomena. In many experiments, such as imaging circulation dynamics, the correlation properties of the observed area are not uniformly distributed, therefore the simulated speckle patterns are supposed to have a desired correlation distribution across the area of interest, which demands the simulation to have a few-to-one mapping [1] between the input matrix and the correlation distribution between simulated speckle patterns.
Various simulation methods have been proposed. The Fast Fourier transform of a phase matrix is a common approach [2], but this method can only synthesize independent speckle patterns. Since many applications of laser speckle are based on their correlation properties, simulation of speckle patterns with a predefined temporal correlation is desirable. Duncan and Kirkpatrick introduced the copula method to generate a sequence of speckle patterns, between which the correlation coefficient could be pre-defined [1,3]. With this method one can investigate the speckle properties of a dynamic object by generating a series of correlated speckle patterns [4,5]. But in this method the generation of the correlated speckle pattern depends on a quantile function and direct Fourier transform, both of which spatially mix the contribution of the input values and the effect of a spatially varying temporal correlation cannot be modelled, i.e. the generated speckle patterns have no location-mapping relationship to the input phase matrices. Therefore this method can only generate speckle patterns with specific temporal correlation coefficients as a whole that are not spatially varying. Although speckle simulation based on the optical transfer function has been investigated [6,7], there is no solution for few-to-one mapping of the correlation coefficients between the input data and the simulated speckle patterns.
In this paper we introduce a simple method based on the concept of the coherent imaging principle, to generate speckle images with a pre-defined correlation matrix to mimic situations with complex correlations across speckle images. The relationship between the input correlation coefficients and the correlation of the simulated speckle patterns are deduced and proved to be power 2.

Generating a variable with a specific correlation to another
Supposing Ω 1 (x,y) and Ω 2 (x,y) are two two-dimensional matrices, which are composed of independent uniformly distributed variables on the interval (-π, π), then e − Ω are similar to the expression of unit optical field at a spatial plane with Ω 1 and Ω 2 as the phases (for monochromatic light) at a moment and W is equivalent to the sum of these two fields' amplitudes. When r is an array, this formula generates a sequence of W where the correlation with 1 i e − Ω is equal to r. This is similar to the temporal correlation changes of speckle images. When r is a two dimensional matrix, the correlation between the corresponding elements in W(x,y) and in 1 i e − Ω are directly related to the values in r, therefore r defines the spatial correlation between W and 1 i e − Ω ; when r changes in three dimensions, it determines the temporal-and-spatial correlation between W and 1 i e − Ω . The variation of the phases Ω 1 and Ω 2 with time and distance need not be considered but could be included when required to simulate a specific optical setup.

Coherent imaging theory in speckle pattern generation
In coherent imaging systems the object plane, the image plane and the coherent transfer function of the imaging system have the following relationship [9]: where Img(x,y) is the Fourier transform of the light field of the image, Obj(x,y) is the Fourier transform of the light field from the object and H(x,y) is the coherent transfer function of the imaging system. In a coherent imaging regime, H(x,y) is the optical pupil function, which is usually a circle or rectangle function depending on the aperture shape. Then the intensity distribution at the image plane, I, can be calculated from the inverse Fourier transform of Img(x,y): where F −1 stands for the inverse Fourier transform.

Synthesis of correlated speckle patterns
Supposing the variable W(x,y) in Eq. (1) is the light field in the object plane, then the intensity distribution of the image plane, which will be a speckle pattern, can be calculated by substituting Obj(x,y) in Eq. (3) with W(x,y) in Eq. (1). The speckle pattern can be expressed as: where F stands for Fourier transform. The introduction of the coherent imaging principle ensures the point-to-point match of element locations between the object and the image. Therefore the correlation of the pixels with the same coordinates in the speckle pattern generated from W and the speckle pattern generated from 1 i e − Ω strictly corresponds to the values defined in the matrix r when based on a large number of random samples. Because Eq. (1) is equivalent to the sum of amplitudes of light fields, the synthesized speckle patterns from Eq. (4) are fully developed and the intensity histogram follows a negative exponential distribution [10]. Since simulation is usually performed with software such as Matlab, the variables here are defined by digital matrices.
It needs to be pointed out that H acts as a low pass filter, therefore in the simulation the central region of H that contains non-zero values can be called the clear aperture and it is smaller than the size of matrix F(W) or the size of W. When speckle patterns are synthesized with Eq. (4), the speckle size can be controlled by the ratio between the size of the clear aperture and the size of matrix W, similar to Kirkpatrick's work [2]. But when this ratio is too small, i.e. the speckle size is big compared to the pixel size, the spatial resolution and the signal-to-noise rate (SNR) of the correlation map calculated from the synthesized speckle images decreases, as observed experimentally. However when the size of the clear aperture is equal to that of W, Eq. (4) does not apply, because there is no frequency filtering and the imaging system is equivalent to free-space light propagation. In this case Eq. (4) approximates to Fraunhofer diffraction supposing that the light propagation distance is much larger than the width of the object plane: This equation has been used to generate speckle patterns, for example see the work of Kirkpatrick et al. [1,2]. Since it is a standard Fourier transform, no spatial correlation mapping is preserved, as expected for a non-imaging free space field propagation. However the correlations between speckle images and the corresponding phase matrix remain.

The correlation between the speckle images
To simplify the equations, 1 W rM r M = + − (6) r can be a one-dimensional, two-dimensional or three-dimensional variable to define the correlation between W and M 1 or the speckle patterns generated from W and M 1 in a pure temporal domain, pure spatial domain or in both spatial and temporal domains. Therefore a sequence of speckle patterns with complicated correlation distribution can be simulated and the correlation between the speckle patterns generated from W and M 1 can be calculated from the following equation: where I W is the intensity of the speckle pattern from W and I M1 is the intensity of the speckle pattern from M 1 . < > denotes the spatial average and σ means the standard deviation.
Because the speckle pattern is fully developed, the PDF of the intensity follows a negative exponential distribution, and the mean is equal to the standard deviation. In the simulation, the total input intensity is kept constant and no energy loss was considered, therefore <I W > = <I M1 > = σ(I W ) = σ(I M1 ). We define this mean as <I>, then Eq. (7) can be reduced to: According to Eq. (1) and Eq. (4), the intensity of the speckle pattern from W can be expressed as: where A M1 and A M2 are the amplitudes of the speckle patterns from M 1 and M 2 , and the superscript * denotes the complex conjugate. Then by substituting Eq. (9) into Eq. (8), we get the correlation: Since the autocorrelation of I M1 is equal to 1, because I M1 and I M2 are independent. In addition, the real part of A M1 and A M2 follows a Gaussian distribution according to the Central Limit Theorem [11] and the mean is equal to zero, therefore the multiplication of the amplitudes follows "product-normal" distribution [12] and the mean remains as 0. Therefore (11) When a non-imaging system was considered, the speckle pattern was synthesized with Eq. (5). The deduction of the relationship between ρ and r was similar to that shown in Eqs. (7)- (11). Using Eq. (5) and Eq. (8), the correlation coefficient ρ can be expressed as: According to the linearity property and cross-correlation theorem of Fourier transform, because M 1 and M 2 are independent, ( ) ( ) ( ) ( ) Therefore Eq. (12) becomes: 2 r ρ = (13) This, together with Eq. (11) proves that the correlation between the intensity of speckle patterns is the square of the input correlation coefficient no matter whether the speckle pattern is synthesized with Eq. (4) or Eq. (5), i.e. we can simulate in the imaging domain.

Simulation results
The simulations were performed with Matlab with the original phase matrices Ω 1 and Ω 2 generated using the function 'RAND' and the values following a uniform distribution on the range (-π, π). In each section below an initial speckle image, I 0 , was synthesized with W = 1 i e − Ω according Eq. (4) or Eq. (5) depending on whether it was simulating free-space propagation or an imaging system for the correlation calculation.

Uniformly correlated speckle patterns in free-space propagation
We defined r = 0, 0.05, 0.1,...,1 and two matrices (600 × 600 pixels) for Ω 1 and Ω 2 . Then the initial image and a sequence of 21 further speckle patterns I 1 , I 2 …I 21 were generated with Eq. (1) and Eq. (5). The correlation coefficients between the sequence of speckle patterns and the initial image were calculated with the following formula: where I 0 is the initial speckle pattern and ρ(0, j) is the correlation coefficients between the initial image and I j , j = 1,2,3,…,21. Cov stands for covariance and σ is the standard deviation. The relationship between ρ and r is plotted in Fig. 1. The calculated correlation coefficients fit well with the curve of r 2 .

Uniformly correlated speckle patterns in coherent imaging system
A matrix H (600 × 600 pixels) was defined as the pupil plane in the frequency domain. The central circular area of H with radius equal to 100 pixels was assigned the value of 1 and the other elements were set to 0, therefore the synthesized speckle size approximated to three image pixels. 21 speckle patterns were generated by using Eq. (1) and Eq. (4) with the same r and Ω 1 , Ω 2 values as in Section 3.1. The correlation coefficients between the initial speckle image and the 21 frames were calculated with Eq. (14). The correlation coefficient as a function of r is shown in Fig. 2(a) and the probability density function (PDF) of the speckle pattern when r was equal to 0.3 is shown in Fig. 2(b). As expected, the square relationship between r and the correlation coefficient of the simulated data remains when the principle of coherent imaging is applied. The PDF of the speckle matches with a negative exponential function as shown in Fig. 2(b), indicating that the speckle patterns are fully developed.

Speckle patterns with complicated correlation distribution and Brownian motion
In this section, r is a matrix containing the desired spatially varying correlation coefficients and the speckle pattern was generated using Eq. (1) and Eq. (4). The pupil function was the same as that in Section 3.1 and 3.2. Figure 3(a) shows the correlation matrix r (600 × 600 pixels) used, which had a value of 0.6 in the 150th to 450th rows and a value of 1 for the other rows. Because the correlation coefficient can only be calculated as a statistical average, we calculated the value of ρ between the speckle pattern generated from W and M 1 by averaging data within the rows. The result is shown in Fig. 3(b), in which the expected values of 0.36 and 1 are marked with a dotted line and a dashed line respectively. Although the curve is noisy, the mean value and the standard deviation are 0.998 and 0.16 respectively in the area with r equal to 1 and 0.347 and 0.096 respectively in the area with r equal to 0.6; i.e. the calculated correlation coefficients are around r 2 as expected. Figure 3(c) shows a demonstration of the correlation distribution. To create this image, a kernel of 11 × 11 pixels was defined and the correlation coefficients were calculated from the pixels in the kernel while the kernel passed through the whole image. This process is similar to the calculation of spatial contrast of speckle images and was accomplished in a Matlab program. The line profile of the 300th column is shown on the left side of Fig. 3(c) to illustrate the correlation values, with levels of 0.36 and 1 indicated. Both Fig. 3(b) and Fig. 3(c) show that the correlation coefficients fluctuate around the estimated values because of the speckle intensity variations and the limited sampling number during the correlation calculation. Further discussion of the signal to noise ratio can be found in the Discussion.  Fig. 3. (a) The phase matrix r used in the simulation, (b) the correlation coefficients of the rows between the speckle images generated from W and M1 (c) the correlation distribution of the two simulated speckle images.
Brownian motion and laminar flow are two common motion types. We synthesized a sequence of speckle pattern frames and the correlation between the first frame and the following ones follows a negative exponential function that corresponds to Brownian motion [13]. The negative exponential correlation is expressed as: where τ c is decorrelation time and τ is the delay time. The 600 × 600 pixel correlation matrix in this simulation is designed as shown in Fig. 4(a) to simulate Brownian motion in the area between the 150th and 450th rows. According to Eq. (13) and Eq. (15), the input correlation values of the pixels between the 150th to 450th rows are defined as: τ c was chosen as 370 μs according to reference [14] and τ was varied from 0 μs to 1.85 ms with an increment of 37 μs. The correlation coefficients were calculated for pixels in the stationary part and in the dynamic part respectively between I 0 and the speckle pattern sequence using Eq. (14). The result is shown in Fig. 4(b) where the asterisks mark the correlation coefficients of the stationary region, the diamonds show the correlation of the Brownian motion region and the line is the negative exponential function. It is clear that the correlation of the stationary parts remains constant while the correlation of the dynamic part decreases as a negative exponential function.

A sequence of speckle patterns with complex correlation distribution
In this section, we simulated speckle patterns with a more complex correlation distribution such as those expected when imaging blood vessels. Note that two approaches were used, the first of which only used these images to demonstrate the feasibility of the simulation method for complex correlation distributions and did not measure or investigate the temporal decorrelation for different areas of the real tissue.
The second approach demonstrated the possibility of creating a simulation based on the correlation features extracted from experimental data. In the first approach, the correlation matrix r 0 was based on an image of a Sprague Dawley rat ear acquired with white light illumination segmented into only five gray-scale levels as shown in Fig. 5(b). We assumed that the larger vessels contained more blood, which strongly absorbed light and therefore exhibited lower intensity in the gray-scale image. It was also expected that these areas would have different decorrelation properties due to the higher flow speed. Therefore five correlation curves could be generated for the five pixel intensity classes by using a theoretically or experimentally derived correlation coefficient. In this simulation the equation Zakharov's work [15] was used, where ρ is the correlation coefficient, x are the acquisition times of the sequence of speckle images, c characterizes the ratio of static part to the total detected light intensity according to reference [15] and it is the lowest correlation coefficient, b is the decorrelation time.
In this simulation we arbitrarily chose c = 0.6, 0.4, 0.3, 0.2, 0.1 and b = 1.5, 1.7, 1.9, 2, 2.2 to assign to the pixels in Fig. 5(b) that correspond to the indicated classes 1 to 5. These values were based on the above assumptions about absorption and flow speed and allowed clearly separated correlation curves and correlation coefficient values to be generated that are reasonable compared to real experiments. Assigning the correlation coefficients for each class in the image generated a correlation matrix for the frame sequence and the correlation curves are shown in Fig. 5(a). Then the independent matrices Ω 1 and Ω 2 were defined to be the same size as the gray-scale image, although this induced poor simulation results due to the statistical residuals, as discussed in the following section. Equation (4) was used to generate the speckle patterns in an imaging domain and 20 frames were synthesized. Figure 5(c) shows the square root of the calculated correlation matrix between the first and the 20th synthesized speckle images. It is noticeable that the intensity distributions of Fig. 5(b) and Fig. 5(c) are similar although Fig. 5(c) exhibits statistical noise. We also calculated the temporal contrast image from the 20 simulated speckle images and the result is shown in Fig. 5(d). The variation in contrast values corresponds well to the chosen correlation coefficients. A second approach was used to further demonstrate the simulation of blood circulation. We experimentally recorded 20 speckle images of a 5 × 3 mm area of a Sprague Dawley rat ear using a 671 nm laser with the camera frame rate equal to 40 fps and the exposure time equal to 20 ms. Then the temporal contrast of the observed area was calculated from the 20 speckle frames and the contrast image is shown in Fig. 6(a). This image was denoised, and segmented into four parts, shown in Fig. 6(b), according to the contrast values. Then the correlation coefficients were calculated respectively within the four contrast classes between the first and the 20 successive speckle frames and are shown in Fig. 6(d). The correlation curve of region 4 within the dashed square in Fig. 6(a) corresponds to the contrast background.
Then a sequence of 20 speckle frames was synthesized, again using the independent phase matrices Ω 1 and Ω 2 and by assigning the correlation between the frames using the data in Fig. 6(d) applied to the regions defined in Fig. 6(b). The simulation was run 10 times and 10 contrast images were calculated and averaged to remove the noise. The result is shown in Fig. 6(c). The repetition of the simulation process was required because the correlation curves were noisy, being based on experimental data rather than purely simulated as in Fig. 5. A comparison of the contrast profile between the experimental data and the simulated data is shown in Fig. 6(e) which shows that the contrast distribution agrees with the segment layout and with the experimental contrast image. The contrast profile of the simulated and the experimental data fit well although the absolute values are different because the experimental images are not fully developed. The normalized error is 5% in the experimental data but is about 40% in the simulated result, giving a lower SNR. However, the normalized error of the simulation can be reduced by decreasing the speckle size and increasing the number of simulations to be averaged. For instance the normalized error was reduced by 50% when the speckle size was decreased to be 2 pixels, or when the number of simulation runs was increased to 30. The boarders of the vessels in Fig. 6(c) are not the same as in Fig. 6(a) because the contrast difference of the simulated result in region 2 and region 3 is not clear due to the low SNR and this can be improved with the methods mentioned above and with less noisy correlation curves. More classes could be used depending on the quality of the raw experimental speckle patterns.

Discussion
Laser speckle analysis has been applied to investigate the changes in blood circulation induced by disease, drugs, and therapeutic response. Although point detection of laser speckle can measure the changes of blood circulation, the investigation of laser speckle properties in an imaging domain can offer additional information. Therefore simulating both temporally and spatially correlated speckle patterns is more desirable than the simulation of speckle patterns that are only temporally correlated. In addition, animal trials are limited by the availability of animal samples and the strict experimental environment. Moreover unproven experimental and data processing methods may induce unnecessary animal sacrifices and loss of time. However the previous simulation methods are not able to simulate this type of speckle pattern. Therefore we proposed this new method to generate speckle patterns with both pre-defined temporal and spatial correlation distribution, which may be experimentally or theoretically derived. With this method laser speckle patterns with complex flow structures can be synthesized.
In this method the square relationship between the input correlation coefficient and the calculated correlation from the simulated speckle patterns is deduced based on the negative exponential distribution property of the speckle intensity, therefore generating fully developed speckle patterns is essential since the intensity of partly developed speckle pattern follow other distributions. There are two conditions to ensure fully developed speckle patterns: the phases Ω 1 and Ω 2 are uniformly distributed on the interval of (-π, π) and the speckle size is no smaller than 2 pixels, which means the size of clear pupil function H(x,y) must be equal to or smaller than half the size of the phase matrices. Sometimes Gaussian distributed variables are used to generate correlated variables. In this case the Gaussian distributed phases need to have a mean equal to 0 and a standard deviation larger than 2π to guarantee the generation of fully developed speckle patterns [10].
The correlation coefficient is a statistical parameter, therefore correlation fluctuationalso called noise -is inevitable in the correlation map calculated from the synthesized speckle patterns. This is due to the small number of spatial sampling points, although it does not affect the feature matching of the correlation distribution between the input data and the simulated data and can be decreased by averaging the result over multiple simulation runs.
When the correlation coefficient distribution was calculated with the kernel method, as in Fig. 3, the noise comes from the statistical error with the minimum determined by the kernel size. This is demonstrated in Fig. 7 with a kernel size equal to 11 × 11.   Fig. 3(a) and the square root of the simulated correlation distribution in Fig. 3(c), which shows the noise distribution. The noise level is lower in the area with r equal to 1, which is reasonable because the pixel intensities in the speckle patterns generated from W and M 1 are identical in the area with r equal to 1 therefore the limitation of sampling number induces less statistical errors than it does in the area with r equal to 0.6. The noise in the area with r equal to 1 comes from the different way of calculating the covariance and standard deviation and it can be removed by programing the calculation of covariance using the same kernel scanning method with the Matlab built-in stdfilt if stdfilt used to calculate the standard deviation. But here we just keep this result since it doesn't change the trend of the noise as a function of r and the number of averaged simulations. Figure 7(b) shows that the standard deviation is lowest when r is equal to 1 and reaches a maximum when r is equal to 0.5, for the same reasons just mentioned. The standard deviation drops slightly as r decreases from 0.5 to zero because the kernel calculation is less affected by the local intensity divergence until the two speckle patterns are totally independent. Figure 7(c) proves that the error can be reduced by averaging over multiple simulations. After 30 runs, the standard deviation is only 0.045, which may be reduced further by increasing the kernel size.
Instead of using a kernel to calculate the correlation, another option is to run multiple simulations with different Ω 1 and Ω 2 to generate multiple pairs of W and Ω 1 speckle images. The correlation can then be calculated for each pixel in the image across the multiple runs, resulting in a higher spatial resolution but requiring more simulation runs to find the correlation. Here we call this method single-pixel based method and is demonstrated in Fig. 8.
We used 200 random Ω 1 and Ω 2 matrices to generate 200 frames respectively from M 1 and W with r the same as the one shown in Fig. 3(a) and calculated the correlation coefficient at each pixel across the 200 frames, as shown in Fig. 8(a). It is clear that Fig.  8(a) has higher spatial resolution than kernel method and the SNR was also tested by subtracting Fig. 8(a) from the input ground truth -the square of Fig. 3(a) -and the result is shown in Fig. 8(b). The regions with r equal to 1 have zero error in this case, since each pair of W and M 1 are identical. The noise is further illustrated as line profiles of the 300th common in Fig. 8(c) for different values of r ( 0.1, 0.2, 0.3, 0.9 ). The noise is lower when r is higher, the same as was found above for the kernel method, and the biggest fluctuation is around ± 0.05. We randomly tested the smallest correlation difference and found that a correlation difference between 0.01 and 0.03 could be resolved in image domain as shown in Fig. 8(d). This limit is expected to reduce with the value of r closer to 1. When using a 2D correlation image extracted from the experimental data as the input correlation matrix and the image is noisy, such as when the values of the surrounding area vary in the same range as those of the features to be detected (e.g. a number of small vessels spreading in the observing area when using a vessel image as the input correlation matrix), the calculated correlation image from the synthesized speckle patterns will lose spatial resolution and the small structures become unperceivable. A demonstration is shown in Fig. 9. In the second simulation in Section 3.4 the image was segmented to allow the correlation coefficients to be more accurately calculated for each class. Ideally the simulated correlation coefficient map would be calculated based on the temporal intensity changes of every single pixel in the speckle image, but this demands a large number of speckle frames. A good camera may also increase the SNR of the correlation coefficient map. This proposed simulation method is therefore able to generate speckle patterns with a specified correlation distribution in the spatial and temporal domains, and also to produce the contrast images. The performance depends on the accuracy with which the correlation map may be derived theoretically or from an experimental result.

Conclusion
In this paper we proposed a new method to simulate a single speckle image or a sequence of speckle images with either uniform or arbitrary correlation distribution in the spatial and temporal domains. The few-to-one mapping of correlation coefficients from the object to the image was demonstrated both by mathematical deduction and simulation. It was proved that the correlation coefficients between the synthesized speckle images are the square of those between the input correlation factors, which enables a simple conversion to any desired correlation mode. This work solved the issue of few-to-one mapping of speckle pattern simulation and provided an effective way to create simulations from in vivo experiments to evaluate processing methods.