Relation between the contrast in time integrated dynamic speckle patterns and the power spectral density of their temporal intensity fluctuations

Matthijs J. Draijer; Erwin Hondebrink; Marcus Larsson; Ton G. van Leeuwen; Wiendelt Steenbergen

doi:10.1364/OE.18.021883

1. Introduction

The speckle phenomenon is widely used for determining tissue perfusion maps [1, 2]. In general, the tissue is illuminated with coherent laser light. A fraction of the laser light interacts with moving red blood cells and obtains a Doppler shift. Doppler shifted and unshifted light which are diffusely scattered from the medium create a dynamic speckle pattern on a plane of observation. Various methods of analyzing this speckle pattern have led to two separate modalities of flux imaging. In the first modality [3] the power spectrum P(ν) of intensity fluctuations generated in the dynamic speckle pattern is analyzed in terms of its moments given by :

M_{i} \equiv \int_{- \infty}^{\infty} ν^{i} P (ν) d ν

where the zeroth order moment (i = 0) is a measure for the concentration of red blood cells and the first order moment (i = 1) is a measure for the flux or perfusion [3]. The physics behind this modality is well-known and it has been shown by Bonner and Nossal [3] that, for low blood concentrations, the concentration of red blood cells and their average velocity are both linearly represented by the power spectral moments of equation 1.

In the second modality, referred to as laser speckle contrast methods, comprising Laser Speckle Contrast Analysis (LASCA) [4] and Laser Speckle Imaging (LSI) [5] the contrast in the speckle pattern is used as a measure for perfusion [2]. In these techniques the changing speckle pattern is averaged over a time interval in the order of the speckle decorrelation-time (i.e. in the millisecond range) leading to speckle blurring. The level of blurring is related to the movement within the illuminated medium and is quantified by the speckle contrast C, which is usually defined as the ratio of the standard deviation σ of the intensity I of the blurred image to the mean intensity 〈I〉 of the speckle pattern :

C \equiv \frac{σ}{〈 I 〉} = \frac{\sqrt{〈 I^{2} 〉 - {〈 I 〉}^{2}}}{〈 I 〉}

where the brackets denote spatial averaging. However, for ergodic speckle patterns, spatial and temporal averaging will give identical results.

This is the case for completely dynamic speckle patterns, hence without static component. Retrieving flux from contrast of time averaged speckle patterns lacks a generally accepted theoretical framework.

As pointed out by Boas and Dunn in their recent review [6], as yet speckle contrast flowmetry modeling has focused on retrieving velocity information of scattering particles rather than their flux which includes concentration [7–9]. The presence of a static speckle component is usually not taken into account, with two recent exceptions [9, 10]. Already Bonner and Nossal [3] take into account the presence of a static tissue matrix, and the experimental evidence given e.g. by Nilsson et al. [11] show that laser Doppler can give a relative measurement of red cell flux within a static environment, at least for a certain range of concentrations.

Furthermore, although less problematic as long as only relative velocities are measured, models for speckle contrast flowmetry are inspired by Dynamic Light Scattering theories, implying assumptions regarding the dynamics of the particles and the associated optical intensity correlations, which is known to lead to model-dependent velocity estimations [8, 12]. Finally they often assume single scattering by moving particles.

Here we present a theory which connects the contrast in time integrated dynamic speckle patterns and the power spectral density of temporal intensity fluctuations of non-integrated speckle patterns without prior assumptions regarding the speed distribution of particles and the extent of multiple scattering. The theory includes speckle patterns with an arbitrary large static component caused by the presence of non-Doppler shifted light. We show that a contrast based parameter which is closest to the weighed spectral moments as proposed by Bonner and Nossal is 1 – C². In this definition, the flux parameter will increase with increasing flow or concentration. Verification of the theory will be done by simulated speckle patterns. Linking speckle contrast flowmetry to a model based on the power spectrum of intensity fluctuations may enable quantification in terms of flux rather than velocity only.

2. Theory

We introduce time dependent function f(t) which is the intensity fluctuation of a pixel, and its moving average $u (t) = \frac{1}{T} \int_{t - T}^{t} f (τ) d τ$ with T the integration time. F(ω) and U(ω) are the Fourier transforms of f(t) and u(t) respectively. From general Fourier transform properties for time domain shifts and integrations [13], U and F are related as

U (T, ω) = \frac{1}{i ω T} [1 - exp (- i ω T)] F (ω)

So the relation between F(ω) and U(ω) is given by the amplitude transfer function

H (T, ω) = \frac{1}{i ω T} [1 - exp (- i ω T)]

with gain

| H (T, ν) | = \frac{1}{T} \sqrt{\frac{1 - cos (2 π ν T)}{2 π^{2} ν^{2}}}

and ν given in Hz. Figure 1 shows an example of |H(T,ν)| as a function of ν for integration times T of 1, 5 and 10 ms. Assume a blurred speckle pattern I(x,t) in which blurring is realized by averaging a dynamic speckle pattern within window T. Assuming statistical homogeneity in time and space, this speckle pattern can be decomposed in an average value 〈I〉, a static, time independent spatial fluctuation I_s(x) and a time- and space dependent fluctuation I_T (x, t) which depends on the integration time T. These components are shown in Fig. 2 which features an example of a blurred intensity I along a line in x-direction. Hence I(x,t) = 〈I〉 + I_s(x) + I_T (x, t). I_s will be nonzero when part of the light producing the speckle pattern has a constant phase, for instance due to interaction with static objects only. For a static-dynamic speckle pattern, for T → ∞ we obtain I_T → 0. For a fully dynamic speckle pattern, I_s(x) = 0 while for T → ∞ also I_T (x, t) → 0 and I(x,t) → 〈I〉. Substituting the above form of I(x,t) in equation 2 gives :

C^{2} = \frac{〈 {(I_{s} + I_{T})}^{2} 〉}{{〈 I 〉}^{2}} = \frac{(〈 I_{s}^{2} 〉 + 〈 I_{T}^{2} 〉)}{{〈 I 〉}^{2}}

where we assumed 〈I_sI_T〉 = 0, since for each possible value of I_s, positive and negative values of I_T occur.

Fig. 1: Gain of transfer function |H(T, ν)| for moving average operation with integration times T = 1, 5 and 10 ms.

Download Full Size | PDF

Fig. 2: Decomposition of the blurred intensity I (solid) into average value 〈I〉 (dotted), static intensity fluctuation I_s (dashed) and time-dependent blurred fluctuation I_T.

Download Full Size | PDF

Furthermore, contrast as a function of integration time T can be written as :

C^{2} (T) = C^{2} (0) + \int_{0}^{T} \frac{d C^{2}}{d \tilde{T}} d \tilde{T}

The intensity in an unblurred polarized speckle pattern will have an exponential probability density function [14], so C²(0) equals unity. Since in the right hand side of equation 3 only I_T (x, t) depends on T, substitution of equation 3 in equation 4 gives

C^{2} (T) = 1 + \frac{1}{{〈 I 〉}^{2}} \int_{0}^{T} \frac{\partial}{\partial \tilde{T}} 〈 I_{T}^{2} 〉 d \tilde{T}

By using Parsevals’s theorem and the transfer function H(T,ν),

〈 I_{T}^{2} 〉

can be written as :

〈 I_{T}^{2} 〉 = \int_{- \infty}^{\infty} P (ν) {| H (T, ν) |}^{2} d ν

Strictly speaking Parsevals theorem implies averaging over all time in the lefthand side. Since static-dynamic speckle patterns are non-ergodic, this should be done both for regions within the speckle pattern that are bright or dark. However, averaging over all space, as indicated by 〈〉, is equivalent to temporal averaging in a limited part of space. Substituting Eq. (6) into Eq. (5) results in :

C^{2} (T) = 1 + \frac{1}{{〈 I 〉}^{2}} \int_{- \infty}^{\infty} P (ν) {[{| H (T, ν) |}^{2}]}_{0}^{T} d ν

which, using Eq. (1) with i = 0, reduces to

C^{2} (T) = 1 - \frac{M_{0}}{{〈 I 〉}^{2}} + \frac{1}{{〈 I 〉}^{2}} \int_{- \infty}^{\infty} P (ν) {| H (T, ν) |}^{2} d ν

since |H(0, ν)| = 1. Equation (8) expresses the speckle contrast in terms of the power spectrum of the local temporal intensity fluctuations in the speckle pattern. For T ↓ 0, it holds that H(T, ν) → 1 for all frequencies, reducing Eq. (8) to C² = 1, which is the required value for a snapshot of the speckle pattern. For T → ∞, H(T, ν) approaches zero and the last term on the right-hand side cancels out, resulting in C²(T → ∞) = 1 – M₀ /〈I〉². For a completely dynamic speckle pattern, the property of ergodicity leads to M₀/〈I〉² = 1 so C²(T → ∞) → 0. Hence, Eq. (8) shows the required behavior for extreme cases.

Since M₀/〈I〉² = f_D(2 – f_D) with f_D the fraction of Doppler shifted light [15], the limiting behavior of Eq. (8) agrees with the form given by Zakharov et al. [10] for quantitative modeling of laser speckle imaging.

To validate Eq. (8), completely dynamic artificial speckle patterns are generated by making use of the concept of a copula [16]. For the speed of change of the speckle pattern, a time scale was chosen which was realistic for tissue speckle. The dynamic speckle pattern was recorded for a total duration of 12.8 seconds.

From the dynamic speckle pattern, for 10 random pixels the intensity as a function of time was extracted. Per pixel, from the time signal the temporal contrast was determined from its definition in Eq. (2). Furthermore in each pixel, from the power spectrum of the time trace the contrast was predicted for different values of T using Eq. (8).

The spatial contrast is determined in concentric regions of 7 × 7 pixels around 10 random pixels. In each region, the contrast was determined from the definition in Eq. (8) as well as predicted by Eq. (2) based on the averaged power spectrum in the region. An example of an averaged power spectrum over 7 × 7 pixels is shown in Fig. 3. The average contrast values and their standard deviations are shown in Fig. 4.

Fig. 3: Power spectrum averaged over an area of 7 × 7 pixels of the artificial speckle pattern used to determine the spatial contrast.

Download Full Size | PDF

Fig. 4: Measured (open symbols) and predicted (closed symbols) speckle contrast values for temporal (squares) and spatial (circles) contrast as a function of integration time T, for artificial speckles.

Download Full Size | PDF

Figure 4 shows that Eq. (8) allows for prediction of both temporal and spatial speckle contrast values based on the power spectrum of the associated intensity fluctuations. There is a much better agreement between the simulated and predicted contrast-curves for the case of temporal contrast than for spatial contrast. Furthermore the error bars for temporal contrast are smaller than for spatial contrast. The discrepancy for the spatial contrast can be explained from the fact that in the limited region of interest of 7 × 7 pixels, the speckle pattern does not exhibit all intensity variations which are present in the complete speckle pattern. For integration times above 1 ms (i.e., which are normally used in speckle contrast techniques [2]) there is good agreement between the predicted and simulated spatial contrast values. The smaller variation of temporal contrast compared to spatial contrast can be explained from the fact that in each of the 10 randomly chosen pixels, spatial contrast is obtained from averaging over 49 pixels in the surrounding region of interest, while the temporal contrast is obtained from all time points. Note the fact that for integration times > 0.5 second still the contrast did not reach zero, presumably due to low frequencies (e.g., below 1 Hz) which are present in the signal.

Given the suitability of the first order power spectral moment of intensity fluctuations as an estimator of particle flux within a static turbid matrix, as shown by Bonner and Nossal [3] for low particle flux, the expression derived in this paper for the contrast in blurred speckle forms the basis for further study of speckle contrast techniques. Here only a first step will be made. Clearly, Eq. (8) shows that speckle contrast provides an integral over the power spectrum P(ν) weighed with |H(T, ν)|². For the contrast based flux parameter 1 – C² we can derive from Eq. (8) that :

1 - C^{2} = \frac{1}{{〈 I 〉}^{2}} \int_{- \infty}^{\infty} (1 - {| H (T, ν) |}^{2}) P (ν) d ν

In Fig. 5 the spectral weighting function in Eq. (9) is shown for integration times of 1, 2, 5 and 15 ms, respectively. The weighting function increases nonlinearly from 0 to 1 at ν = 1/T. For ν > 1/T the weighting function is almost constant, showing decaying oscillations between 1 and 0.95. Hence, for integration times 5 ms < T < 15 ms realistic for speckle contrast techniques, frequency dependent weighting is only performed in a frequency interval between 0 Hz and 67 Hz < ν < 200 Hz. For higher frequencies, 1 – C² mainly provides the zero order moment of the power spectrum. For realistic integration times and for low concentrations of moving particles that are assumed in the theory of Bonner and Nossal [3], speckle contrast mainly provides information regarding the concentration of moving particles rather than their speed. Information about speed variations is only conveyed inasmuch as these variations affect the spectral broadening in the interval 0 < ν < 1/T (in Hz). In this interval the frequency weighting is nonlinear, with an approximately 2^nd order weighting (∝ ν²) for ν → 0. The second order weighting can be extended to higher frequencies by reducing integration time T. For instance, assuming that the spectral weighting function associated with 1 – C² has a second order behavior for 0 < ν < 1/2T, the power spectrum shown in Fig. 3, with a width of 5 kHz which is typical for physiological perfusion, would be obtained for an integration time of approximately 0.1 ms. This will give a contrast value which is only slightly smaller than 1, and therefore will not be very sensitive. Figure 5 also suggests that speckle contrast methods will be particularly sensitive to changes in the power spectrum in the low frequency range, e.g. caused by overall tissue motion or speed variations of particles moving at a low speed. For other flux estimators which provide higher flux or velocity values for lower contrast, such as 1/C, a similar analysis may be made, however this will be mathematically less elegant.

Fig. 5: Spectral weighting functions realized by 1 – C², with C the contrast of integrated speckle for integration times T of 1, 2, 5 and 15 ms.

Download Full Size | PDF

3. Conclusion

In this paper, we have presented a theory which expresses the contrast in time integrated dynamic speckle patterns in terms of the power spectral density of their local temporal intensity fluctuations. The theory covers mixed static-dynamic speckle patterns, provided that they are statistically homogeneous. Verification was done on computer-created fully dynamic speckle patterns.

We show that with speckle contrast C, the flux parameter 1 – C² provides a weighted average of the power spectrum of photocurrent fluctuations similar to that in laser Doppler flowmetry, however with nonlinear rather than linear weighting. For very small integration times 1 – C² would imply a second order weighting. For realistic integration times T > 5 ms, speckle contrast mainly provides information regarding the concentration of particles moving within a static matrix, with speed information only present as far as represented by the power spectrum for frequencies between zero and 1/T Hz. The presented theory will enable further research into the use of speckle contrast as an estimator of tissue perfusion.

References and links

1. P. Vennemann, R. Lindken, and J. Westerweel, “In vivo whole-field blood velocity measurement techniques,” Exp. Fluids 42(4), 495–511 (2007). [CrossRef]

2. M. J. Draijer, E. Hondebrink, T. G. van Leeuwen, and W. Steenbergen, “Review of laser speckle contrast techniques for visualizing tissue perfusion,” Lasers Med. Sci. 24(4), 639–651 (2009). [CrossRef]

3. R. Bonner and R. Nossal, “Model for Laser Doppler Measurements of Blood Flow in Tissue,” Appl. Opt. 20(12), 2097–2107 (1981). [CrossRef] [PubMed]

4. J. D. Briers and S. Webster, “Laser speckle contrast analysis LASCA): a nonscanning, full-field technique for monitoring capillary blood flow,” J. Biomed. Opt. 1, 174–179 (1996). [CrossRef]

5. H. Cheng, Q. Luo, S. Zeng, S. Chen, J. Cen, and H. Gong, “Modified laser speckle imaging method with improved spatial resolution,” J. Biomed. Opt. 8(3), 559–564 (2003). [CrossRef] [PubMed]

6. D. A. Boas and A. K. Dunn, “Laser speckle contrast imaging in biomedical optics,” J. Biomed. Opt. 15(1), 011109 (2010). [CrossRef] [PubMed]

7. J. D. Briers and S. Webster, “Quasi real-time digital version of single-exposure speckle photography for full-field monitoring of velocity or flow fields,” Opt. Commun. 116, 36–42 (1995). [CrossRef]

8. J. C. Ramirez-San-Juan, R. Ramos-Garcia, I. Guizar-Iturbide, G. Martinez-Niconoff, and B. Choi, “Impact of velocity distribution assumption on simplified laser speckle imaging equation,” Opt. Express 16(5), 3197–3203 (2008). [CrossRef] [PubMed]

9. A. B. Parthasarathy, W. J. Tom, A. Gopal, X. Zhang, and A. K. Dunn, “Robust flow measurement with multi-exposure speckle imaging,” Opt. Express 16(3), 1975–1989 (2008). [CrossRef] [PubMed]

10. 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]

11. G. E. Nilsson, T. Tenland, and P. A. Öberg, “A New Instrument for Continuous Measurement of Tissue Blood FlowBy LightBeating Spectroscopy,” IEEE Trans. Biomed. Engin. 27(1), 12–18 (1980). [CrossRef]

12. 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]

13. P. A. Lynn, Electronic signals and systems (Macmillan, London, 1986).

14. J. W. Goodman and G. Parry, Laser Speckle and Related Phenomena (Springer-Verlag, New York, 1975).

15. A. Serov, W. Steenbergen, and F. F. M. de Mul, “Prediction of the photodetector signal generated by Doppler-induced speckle fluctuations: theory and some validations,” J. Opt. Soc. Am. A 18(3), 622–630 (2001). [CrossRef]

16. D. D. Duncan, S. J. Kirkpatrick, and R. K. Wang, “Statistics of local speckle contrast,” J. Opt. Soc. Am. A Opt. Image Sci. Vis. 25(1), 9–15 (2008). [CrossRef]

Relation between the contrast in time integrated dynamic speckle patterns and the power spectral density of their temporal intensity fluctuations

Abstract

1. Introduction

2. Theory

3. Conclusion

References and links

Cited By

Figures (5)

Equations (12)

Optics Express