Simulation of speckle patterns with pre-defined correlation distributions

Lipei Song; Zhen Zhou; Xueyan Wang; Xing Zhao; Daniel S. Elson

doi:10.1364/BOE.7.000798

1. 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.

2. Theory

2.1 Generating a variable with a specific correlation to another

Supposing Ω₁(x,y) and Ω₂(x,y) are two two-dimensional matrices, which are composed of independent uniformly distributed variables on the interval (-π, π), then $e^{- i Ω_{1}}$ and $e^{- i Ω_{2}}$ are uncorrelated. One method to generate a new variable W with a specific correlation to $e^{- i Ω_{1}}$ and $e^{- i Ω_{2}}$ is as follows [8]:

W = r e^{- i Ω_{1}} + \sqrt{1 - r^{2}} e^{- i Ω_{2}}

where r is the correlation coefficient factor between W and

e^{- i Ω_{1}}

. When r is adjusted from 1 to 0, W transitions from fully dependent on

e^{- i Ω_{1}}

to fully independent of

e^{- i Ω_{1}}

. It can be noticed that

e^{- i Ω_{1}}

and

e^{- i Ω_{2}}

are similar to the expression of unit optical field at a spatial plane with Ω₁ and Ω₂ 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

e^{- i Ω_{1}}

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

e^{- i Ω_{1}}

are directly related to the values in r, therefore r defines the spatial correlation between W and

e^{- i Ω_{1}}

; when r changes in three dimensions, it determines the temporal-and-spatial correlation between W and

e^{- i Ω_{1}}

. The variation of the phases Ω₁ and Ω₂ with time and distance need not be considered but could be included when required to simulate a specific optical setup.

2.2 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]:

Img (x, y) = O b j (x, y) H (x, y)

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):

I = {| F^{- 1} (Img (x, y)) |}^{2} = {| F^{- 1} [O b j (x, y) H (x, y)] |}^{2}

where F ⁻¹ stands for the inverse Fourier transform.

2.3 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:

I = {| F^{- 1} [F (W) H (x, y)] |}^{2}

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

e^{- i Ω_{1}}

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:

I = {| F (W) |}^{2}

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.

2.4 The correlation between the speckle images

To simplify the equations, $e^{- i Ω_{1}}$ is denoted by M ₁ and $e^{- i Ω_{2}}$ by M ₂. The correlation between W and M ₁ is r and M ₁ is independent of M ₂. Then Eq. (1) becomes:

W = r M_{1} + \sqrt{1 - r^{2}} M_{2}

r can be a one-dimensional, two-dimensional or three-dimensional variable to define the correlation between W and M ₁ or the speckle patterns generated from W and M ₁ 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₁ can be calculated from the following equation:

ρ = \frac{〈 (I_{W} - 〈 I_{W} 〉) (I_{M 1} - 〈 I_{M 1} 〉) 〉}{σ (I_{W}) σ (I_{M 1})}

where I_W is the intensity of the speckle pattern from W and I_M ₁ is the intensity of the speckle pattern from M ₁. < > 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_M ₁> = σ(I_W) = σ(I_M ₁). We define this mean as <I>, then Eq. (7) can be reduced to:

ρ = \frac{〈 I_{W} I_{M 1} 〉 - {〈 I 〉}^{2}}{{〈 I 〉}^{2}}

According to Eq. (1) and Eq. (4), the intensity of the speckle pattern from W can be expressed as:

〈 I_{M 1} (A_{M 1} A_{M 2}^{*} + A_{M 2} A_{M 1}^{*}) 〉 = 0

where A_M ₁ and A_M ₂ are the amplitudes of the speckle patterns from M ₁ and M ₂, and the superscript * denotes the complex conjugate. Then by substituting Eq. (9) into Eq. (8), we get the correlation:

\begin{array}{l} ρ = \frac{〈 [r^{2} I_{M 1}^{} + (1 - r^{2}) I_{M 2}^{} + r \sqrt{1 - r^{2}} (A_{M 1} A_{M 2}^{*} + A_{M 1}^{*} A_{M 2})] I_{M 1} 〉 - {〈 I 〉}^{2}}{{〈 I 〉}^{2}} \\ = \frac{r^{2} 〈 I_{M 1}^{2} 〉 + (1 - r^{2}) 〈 I_{M 1} I_{M 2} 〉 + r \sqrt{1 - r^{2}} 〈 I_{M 1} (A_{M 1} A_{M 2}^{*} + A_{M 1}^{*} A_{M 2}) 〉 - {〈 I 〉}^{2}}{{〈 I 〉}^{2}} \end{array}

Since the autocorrelation of I_M1 is equal to 1, $〈 I_{M 1}^{2} 〉 = 2 {〈 I 〉}^{2}$ according to Eq. (8). Similarly $〈 I_{M 1} I_{M 2} 〉 = {〈 I 〉}^{2}$ because I_M ₁ and I_M ₂ 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 $〈 I_{M 1} (A_{M 1} A_{M 2}^{*} + A_{M 2} A_{M 1}^{*}) 〉 = 0$ and Eq. (10) becomes:

ρ = r^{2}

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:

\begin{array}{l} ρ = \frac{r^{2} 〈 I_{M 1}^{2} 〉 + r \sqrt{1 - r^{2}} 〈 I_{M 1} (F (M_{1}) F^{*} (M_{2}) + F * (M_{1}) F (M_{2})) 〉}{{〈 I 〉}^{2}} \\ + \frac{(1 - r^{2}) 〈 I_{M 1} I_{M 2} 〉}{{〈 I 〉}^{2}} - 1 \end{array}

According to the linearity property and cross-correlation theorem of Fourier transform, because M ₁ and M ₂ are independent, $F (M_{1}) F^{*} (M_{2}) = F^{*} (M_{1}) F (M_{2})$ = 0. Therefore Eq. (12) becomes:

ρ = r^{2}

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.

3. Simulation results

The simulations were performed with Matlab with the original phase matrices Ω₁ and Ω₂ generated using the function ‘RAND’ and the values following a uniform distribution on the range (-π, π). In each section below an initial speckle image, I₀, was synthesized with W = $e^{- i Ω_{1}}$ according Eq. (4) or Eq. (5) depending on whether it was simulating free-space propagation or an imaging system for the correlation calculation.

3.1 Uniformly correlated speckle patterns in free-space propagation

We defined r = 0, 0.05, 0.1,...,1 and two matrices (600 × 600 pixels) for Ω₁ and Ω₂. Then the initial image and a sequence of 21 further speckle patterns I₁, I₂…I₂₁ 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:

ρ (0, j) = \frac{C o v (I_{0}, I_{j})}{σ (I_{0}) σ (I_{j})}

where I₀ 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 ².

Fig. 1 The relationship between correlation coefficient of speckle patterns and that of the phase matrices.

Download Full Size | PDF

3.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 Ω₁, Ω₂ 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.

Fig. 2 (a) The correlation coefficient of the speckle patterns as a function of r; (b) the probability density function of the speckle pattern.

Download Full Size | PDF

3.3 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 ₁ 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 ² 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 M₁ (c) the correlation distribution of the two simulated speckle images.

Download Full Size | PDF

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:

ρ = \exp (- τ / τ_{c})

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:

r = \sqrt{\exp (- τ / τ_{c})}

τ_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₀ 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.

Fig. 4 (a) The flow model. (b) The simulated correlation for the two regions as a function of the frame index. The asterisks mark the correlation coefficients of the stationary area with input values equal to 1; the diamonds represent the correlation coefficients in the Brownian motion region, and the line is a negative exponential function.

Download Full Size | PDF

3.4 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₀ 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 $ρ = (1 - c) \exp (- x / b) + c$ derived from P. 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.

Fig. 5 (a) The correlation coefficients of different segments and different frame indices; (b) the correlation distribution defined by five gray levels. This correlation image and the correlation curves were used to synthesize the 20 speckle patterns; (c) the correlation image calculated from the 1st and the 20th simulated speckle patterns; (d) temporal contrast calculated from the 20 simulated speckle patterns.

Download Full Size | PDF

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 Ω₁ and Ω₂ 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.

Fig. 6 (a) Denoised temporal contrast from the experiment; (b) four labelled regions based on the temporal contrast image shown in (a); (c) the temporal contrast image calculated from the simulation; (d) correlation coefficients used in the simulation; (e) comparison of the contrast profile along the pixels marked by the dashed line in (c) from the experiment and the simulation after smoothing.

Download Full Size | PDF

Then a sequence of 20 speckle frames was synthesized, again using the independent phase matrices Ω₁ and Ω₂ 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.

4. 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 Ω₁ and Ω₂ 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 fluctuation - also 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. 7 Investigation of noise level of the correlation calculated with a kernel size equal to 11 × 11. (a) The subtraction of Fig. 3(a) and the square root of Fig. 3(c); (b) the standard deviation of the simulated correlation coefficients as a function of r; (c) the standard deviation of the simulated correlation as a function of the number of correlation images that are averaged.

Download Full Size | PDF

Figure 7(a) is the subtraction of the ground truth data in 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 ₁ 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 Ω₁ and Ω₂ to generate multiple pairs of W and Ω₁ 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 .

Fig. 8 Single-pixel based correlation coefficient calculated from simulated speckle images. (a) the correlation map when using the matrix shown in Fig. 3(a); (b) the subtraction of the square of Fig. 3(a) and Fig. 8(a); (c) the line profile of the correlation coefficients along the 300th column with different input value of correlation coefficients; (d) the smallest correlation coefficient difference that can be shown in the image domain. In this case the input values of r are equal to 0.01 and 0.03 respectively for the two areas.

Download Full Size | PDF

We used 200 random Ω₁ and Ω₂ matrices to generate 200 frames respectively from M ₁ 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 ₁ are identical. The noise is further illustrated as line profiles of the 300th common in Fig. 8(c) for different values of r ( $\sqrt{0.1}, \sqrt{0.2}, \sqrt{0.3}, \sqrt{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 .

Fig. 9 Demonstration of a simulation using noisy images as the input correlation matrix. (a) the input correlation matrix; (b) the correlation map calculated from synthesized speckle images; (c) the temporal contrast image calculated from the synthesized speckle images.

Download Full Size | PDF

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.

5. 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.

Acknowledgments

This work was sponsored by the open project of State Key Lab KF201301, National Natural Science Foundation of China (NSFC) (11474169) and Science and Technology Innovation Program of the Chinese Academy of Sciences. Research funding is also gratefully acknowledged from the ERC grant StG 242991 (Dr. Daniel Elson).

References and links

1. D. D. Duncan and S. J. Kirkpatrick, “The copula: a tool for simulating speckle dynamics,” J. Opt. Soc. Am. A 25(1), 231–237 (2008). [CrossRef] [PubMed]

2. S. J. Kirkpatrick, D. D. Duncan, and E. M. Wells-Gray, “Detrimental effects of speckle-pixel size matching in laser speckle contrast imaging,” Opt. Lett. 33(24), 2886–2888 (2008). [CrossRef] [PubMed]

3. D. D. Duncan and S. J. Kirkpatrick, “Algorithms for simulation of speckle (laser and otherwise),” Biomedical Optics (BiOS) (Optical Society of America, 2008), pp 685505.

4. O. B. Thompson and M. K. Andrews, “Tissue perfusion measurements: multiple-exposure laser speckle analysis generates laser Doppler-like spectra,” J. Biomed. Opt. 15(2), 027015 (2010). [CrossRef] [PubMed]

5. H. Zhang, P. Li, N. Feng, J. Qiu, B. Li, W. Luo, and Q. Luo, “Correcting the detrimental effects of nonuniform intensity distribution on fiber-transmitting laser speckle imaging of blood flow,” Opt. Express 20(1), 508–517 (2012). [CrossRef] [PubMed]

6. V. Nascov, C. Samoilă, and D. Ursuţiu, “Fast computation algorithms for speckle pattern simulation,” AIP Conf. Proc. 1564, 217–222 (2013). [CrossRef]

7. H. Fujii, J. Uozumi, and T. Asakura, “Computer simulation study of image speckle patterns with relation to object surface profile,” J. Opt. Soc. Am. 66(11), 1222–1236 (1976). [CrossRef]

8. D. M. Manfred Gilli and E. Schumann, Numerical Methods and Optimization in Finance (Academic Press, 2011).

9. J. Goodman, Introduction to Fourier Optics (Roberts and Company Publishers, 2004).

10. J. W. Goodman, Speckle Phenomena in Optics: Theory and Applications (Roberts & Co, 2007).

11. J. W. Goodman, Statistical Optics, 2nd ed. (Wiley, 2015).

12. 'Normal Product Distribution', Mathworld.

13. D. D. Duncan and S. J. Kirkpatrick, “Can laser speckle flowmetry be made a quantitative tool?” J. Opt. Soc. Am. A 25(8), 2088–2094 (2008). [CrossRef] [PubMed]

14. V. Rajan, B. Varghese, T. G. van Leeuwen, and W. Steenbergen, “Speckle size and decorrelation time; space–time correlation analysis of coherent light dynamically scattered from turbid media,” Opt. Commun. 281(6), 1755–1760 (2008). [CrossRef]

15. P. Zakharov, A. Völker, A. Buck, B. Weber, and F. Scheffold, “Quantitative modeling of laser speckle imaging,” Opt. Lett. 31(23), 3465–3467 (2006). [CrossRef] [PubMed]

Simulation of speckle patterns with pre-defined correlation distributions

Abstract

1. Introduction

2. Theory

2.1 Generating a variable with a specific correlation to another

2.2 Coherent imaging theory in speckle pattern generation

2.3 Synthesis of correlated speckle patterns

2.4 The correlation between the speckle images

3. Simulation results

3.1 Uniformly correlated speckle patterns in free-space propagation

3.2 Uniformly correlated speckle patterns in coherent imaging system

3.3 Speckle patterns with complicated correlation distribution and Brownian motion

3.4 A sequence of speckle patterns with complex correlation distribution

4. Discussion

5. Conclusion

Acknowledgments

References and links

Cited By

Figures (9)

Equations (16)

Biomedical Optics Express