Tissue perfusion measurements: multiple-exposure laser speckle analysis generates laser Doppler–like spectra

Perfusion levels in the skin microvasculature can be measured and mapped by the related noncontact optical methods of laser Doppler^1- 3 and laser speckle contrast measurements^4- 5 in a variety of configurations. These techniques are applicable to various biomedical imaging tasks such as quantifying the progress of peripheral vascular disease in diabetes or monitoring the reperfusion of skin flaps in plastic surgery.

Doppler perfusion measurements have achieved some medical acceptance, but there has been debate over the correct interpretation of laser speckle perfusion measurements. This paper shows that a multiple-exposure speckle method provides exactly equivalent information to a Doppler measurement.

In both laser Doppler and laser speckle contrast analysis, the tissue under test is illuminated by a laser. A dynamic laser speckle pattern is formed by light multiply scattered in the tissue and returned to an electronic sensor.

If the sensor records the intensity of the field at the detection point directly, with no imaging optics, the configuration is referred to as objective speckle⁶ or far-field speckle⁷; if the tissue or target is imaged by a lens, it is the lens aperture that controls the speckle, and the situation is termed subjective speckle⁶ or image speckle.⁷ It can be shown that the speckle statistics of both methods are generally identical.⁷

A small coherent region of an objective speckle pattern is collected by a photodiode in the case of most Doppler methods, and a large area of a subjective speckle pattern is recorded using a digital camera in laser speckle contrast imaging.

In laser Doppler analysis, the fluctuations in intensity within a single speckle are interpreted as being the result of beating between Doppler-shifted light, which has encountered a moving blood cell in its multiply scattered path through the tissue, and non-Doppler-shifted light, which has encountered only static scatterers. Laser Doppler methods generally calculate a perfusion index, defined as the first moment of the power spectrum of these fluctuations.⁸ This perfusion index has been shown by both calculation and experiment to be proportional to both the concentration and mean speed of moving blood cells in tissue.⁹ Doppler methods record a continuous time series, or the related spectrum, at each measurement point. This reduces the applicability of Doppler methods to imaging tasks, as in order to generate an image, the measurement point must be scanned over the imaging area and sufficient time allowed at each point to measure a power spectrum and the related perfusion index. A typical commercial laser Doppler system takes up to $5min$ to record a $256×256pixel$ image.¹⁰

Typical laser speckle perfusion measurements take a single image of the skin surface, illuminated with laser light, and calculate the perfusion in each area of the image as a function of the speckle contrast in that area. Speckle contrast, defined as the ratio of the standard deviation of the intensity to the mean, falls as increasing blood flow generates more rapid fluctuations in the biospeckle pattern and more blurring at a fixed camera exposure. Speckle measurements are ideally suited to imaging tasks. They generate a full-field perfusion image for each exposure, and their history in such applications is long,¹¹ but their uptake required the introduction of modern high-speed digital cameras. Acceptance has also been hindered by confusion and uncertainty about the quantitative interpretation of the images. It has been argued for many years that laser speckle contrast and laser Doppler measurements are manifestations of the same physical phenomenon, with different methods and interpretation.⁷ While this view has some acceptance, there are still distinctions made between the two methods.

We have proposed that multiexposure speckle contrast measurements should recover the same spectral information that is provided by Doppler methods.¹² In this paper, we show that this is the case, using both analysis of computer-synthesized laser speckle data and experimental results directly comparing laser Doppler and laser speckle contrast analysis. The experimental measurements are made on both the Brownian motion in a small vial of milk and on the blood flow in a finger. Subject to the caveat that all experimental measurements of perfusion require an instrument proportionality factor before perfusion can be expressed in absolute units, speckle contrast measurements are no less quantitative than Doppler-based methods.

Multiple-exposure speckle contrast analysis uses the increase in blurring of dynamic speckle with increasing exposure to infer the frequencies present at each point in the image. For a dynamic speckle image at an infinitesimal exposure, or an exposure shorter than the period of the shortest intensity fluctuation, there is no blurring, and the contrast of fully developed speckle is 1. Increasing the exposure to a particular finite value $Tf$ blurs those fluctuations that have periods shorter than $Tf$, while fluctuations longer than $Tf$ remain unchanged. The speckle contrast is reduced from 1 as the speckle pattern blurs, and the degree of reduction in contrast depends on both $Tf$ and the frequency spectrum of the intensity fluctuations. This dependency is developed more rigorously in the following mathematical section.

Goodman¹³ derives the following equation linking the spatial variance $σ2$ of a speckle pattern captured using a camera with exposure time $T$ with $Ct(τ)$, the temporal autocovariance of the intensity at a point:Display Formula

Previous workers, beginning with Fercher and Briers⁴ have assumed a form of the autocorrelation function. Equation 6 shows how the autocorrelation function can be measured using speckle contrast methods. The corresponding power spectrum can then be found using the Wiener–Khintchine theorem.

The nature of dynamic speckle from a volume scatterer like skin, in which diffuse reflection involves multiple scattering and some scattering paths involve moving scatterers, is quite different from the translating speckle pattern produced by moving a simple single scattering surface. The former appears to “boil” as the speckle pattern evolves with time. This can be modeled by allowing the phase of points in a field to evolve randomly and generating the resulting speckle pattern using Fourier optics.

Simulated speckle data were generated following the methods described by Duncan and Kirkpatrick.⁶ Their algorithms for simulating speckle in MATLAB code are available.¹⁴ To generate a single frame of speckle, we start from a matrix of size $L$ by $L$ containing a square of size $L′$ by $L′$ filled with complex numbers with unity amplitude and randomly distributed phase. This set of vectors represents the phases of randomly scattered light that will form the speckle pattern—or the phase changes applied to a coherent source by the effect of scattering in tissue. To produce the speckle pattern, the $L×L$ 2-D Fourier transform of the random phase matrix is taken and then multiplied pointwise by the complex conjugate. The ratio $L∕L′$ must be 2 or higher in order to simulate a speckle sampling scheme that meets the Nyquist sampling criterion.¹⁵ Under these conditions, the speckle patterns show the correct statistics for polarized, fully developed speckle.⁶

A succession of speckle image frames in which the speckle evolved with time was generated. These frames were stacked to produce a 3-D set of data, termed the speckle cube. Both Doppler and multiexposure speckle contrast were used to calculate the power spectral density (PSD) of the light intensity at a typical point in the frame.

Duncan and Kirkpatrick describe using a copula to generate speckle with defined interframe correlation.¹⁶ We use a simpler technique: to generate speckle frames with a small interframe decorrelation, we add a small, normally distributed random phase to each of the vectors representing a scattered light component and calculate a new speckle frame as described earlier using the new scattered light matrix. Repeating this process for the number of frames required produces a volume of speckle data, as shown in Fig. 1. A single speckle frame, as seen on the smaller square face of the speckle volume, shows typical speckle. The sides of the speckle volume plotted in 3-D represent the change of intensity of a typical single line in the image with time.

For the following results, a set of speckle data with 1024 frames, each $100×100pixels$ with $L∕L′=2$, was used. Each frame is strongly correlated with the immediately adjacent frames, with mean $R$ between adjacent frames of 0.96. This computed decorrelation was set by adjusting the magnitude of the random phase addition described earlier. We define the time step as $1ms$ for convenience and to make the simulation comparable to skin measurements.

Doppler perfusion measurements are typically based on a power spectrum of photodetector current.⁸ Taking the square of the absolute value of a fast Fourier transform (FFT) of intensity in the time dimension at a single spatial point in our synthetic data gives us the power spectral density (PSD) of intensity or detector current at that point. The mean of these spectra over all spatial points in the array gives us a mean PSD using all of the synthetic speckle data, reducing noise.

This PSD analysis is a fair approximation to a Doppler system that samples entirely within a single speckle. In order to find the PSD produced by a Doppler system in which the sensor is larger than the minimum speckle size, we can measure the mean PSD of an integrated area several pixels wide. Figure 2 compares the original PSD with a PSD calculated after first taking the mean intensity over $4×4pixel$ areas. The shape of the spectrum is maintained, as is the DC level; taking an area larger than the minimum speckle size simply reduces the AC components of the signal.

The basis of the multiexposure speckle contrast analysis is a $K(T)$ curve—measurements of speckle contrast $K$ at a range of exposures $T$ (Ref. ¹²). We start the analysis of the simulated speckle by generating one of these curves. Longer exposures in a camera integrate the speckle intensity over time, and this is simulated by taking intensity sums of pixels over a number of adjacent frames. It is convenient to use exposures that increase exponentially, and here the exposure length was increased by increasing the number of summed speckle frames by powers of two. A single frame, with $K=1$, is considered to have exposure $T=0$, and $N$ summed frames to have exposure $T=(N−1)ms$. Any inaccuracies produced by quantization errors in this procedure will be most significant at short exposures, corresponding to high frequencies in the eventual PSD plot. Again, all of the data available was used for smoothing, by calculating the mean contrast $K$ for 1024 single frames $(T=0)$, the mean $K$ of 512 sums of two frames $(T=1ms)$, and so on, with the entire frame used as the speckle calculation region. The result, plotted in Fig. 3, resembles the plots produced in the past from skin measurements.¹² There are some differences, resulting from the fact that the simple speckle-generating algorithm does not reflect speckle evolution in real tissue produced by the actual scatter process and scatterer speeds. These differences will generate a different power spectral density than that from perfused tissue, but the comparison between speckle and Doppler techniques for measuring this PSD remains valid.

A curve was fitted to the calculated points using a spline function. From the smooth contrast versus exposure curve, the autocorrelation function was numerically computed using the double differential equation [Eq. 6] earlier. The autocorrelation function produced is plotted in Fig. 4, together with the temporal autocorrelation function from each pixel calculated directly using the MATLAB routine $xcov$ and then averaged over all pixels to reduce noise.

The two curves are essentially the same. Since the temporal PSD at a pixel may be calculated via the Wiener-Khintchine relation from the temporal autocorrelation function, the PSD computed via speckle contrast must now be expected to equal that computed from an FFT of the photocurrent at a point.

This PSD calculated from the speckle-derived autocorrelation is shown in Fig. 5, with the PSD previously calculated in the Doppler analysis for comparison. The power spectra calculated using the two methods are clearly consistent. The main difference is at high frequencies corresponding to short exposures in the contrast versus exposure curve, where the simulation is most susceptible to quantization errors.

The PSD was measured on both the skin of the volar surface of the forefinger and a test target using both the laser Doppler and laser speckle systems. Measurements were made using both systems sequentially. The Brownian motion in homogenized low-fat milk in a $10-mm$-diam transparent plastic tube, with the surface slightly roughened using fine sandpaper to reduce specular reflections, provided a suitably consistent test target.

The simple laser Doppler system used consisted of a photodiode sensor (Centronics AEPX65), a transimpedance amplifier giving $42MV∕A$, and an HP spectrum analyzer (HP 3589A). The laser, a $658-nm$, $50-mW$ thermo-electrically stabilized single-mode laser diode module from WorldStarTech, was focused to a spot with effective diameter approximately $0.5mm$, and the photodiode placed at $250mm$ from the target. The minimum objective speckle size at this range is about $0.4mm$, compared to the $0.85-mm$ diameter of the photodiode. Despite this mismatch, which reduced the recorded AC signal to some degree, there was a sufficient AC signal to record a Doppler spectrum. Doppler spectra were recorded on both the skin and milk dynamic targets and a static, multiply scattering Teflon target as a control for amplifier, dark-current, and background light noise. This small background measurement was subtracted from the dynamic target spectra. The Doppler spectra are plotted later in Figs. 9,12 as discrete points rather than continuous curves, as the values were recorded from the spectrum analyzer using a cursor function.

Speckle contrast was measured using the same laser as used for the Doppler measurements, with the beam expanded using a lens to cover the target, a monochrome digital industrial camera (Sony XCD-SX910), and custom software to process the data. The camera was fitted with a $75-mm$ lens, a bandpass interference filter matched to the laser wavelength, and a polarizing filter. The polarizing filter is used to restore full contrast, as the multiply scattering target produces a contrast reduction of $1∕√2$ due to the presence of orthogonal interference patterns. The polarizing filter also removes any remaining specular reflections from the surface of the tube. The camera aperture was set to f/16, giving a minimum speckle size at the sensor of $15μm$, larger than the Nyquist minimum size for this sensor of $2×4.65μm$ (Ref. ¹⁵). Speckle contrast calculation regions of $50×50pixels$ and a range of camera exposures from $0.1to100ms$ on the milk target and $0.05to400ms$ on skin were used.

Speckle measurements were made on the Brownian motion milk target using both the conventional imaging setup and an objective speckle setup. In the latter configuration, the lens was removed from the camera, the laser was focused to a spot with effective size approximately $0.5mm$, and the objective speckle pattern generated was recorded directly at the CCD chip. The minimum speckle size in this configuration was approximately 10 times the pixel size. The bandpass and polarizing filters remained in the system, in front of the CCD chip. Both of these laser speckle configurations, and the laser Doppler configuration, are illustrated in Fig. 6.

The laser speckle contrast results for the lensed subjective setup on the Brownian motion milk target are shown in Fig. 7. The speckle contrast increases with reducing exposure, with the expected sigmoidal shape when plotted to be logarithmic in time. Interpolating polynomial functions are fitted to this data set and to the equivalent objective speckle data. The latter data are very similar to the subjective speckle measurements and are omitted from the plots here for clarity.

The temporal autocorrelation function for the objective speckle milk measurements, calculated according to the differential Eq. 6 earlier, is plotted in Fig. 8. The troughs adjacent to the central peak are likely artifacts of our numerical processing and might be removed with a larger number of measured points, and a closer approach of the measured points to zero exposure.

The PSD for milk, computed from the autocorrelation functions found in both objective and subjective speckle experiments, are plotted in Fig. 9, together with the Doppler measurements for comparison. These values are arbitrarily scaled in order to overlay each other, since there is an effectively arbitrary gain value in both systems. The Doppler and speckle spectra are equivalent.

The speckle contrast measurements on skin are plotted in Fig. 10, and the corresponding autocorrelation function in Fig. 11. The PSD estimates by Doppler and speckle for skin are plotted in Fig. 12—again, the Doppler and speckle spectra are equivalent. The small difference between the PSD curves at high frequencies corresponds to the region of the speckle contrast curve extrapolated toward zero exposure, and it is possible that this extrapolation generated the small error.

It is clear from the preceding results that the same spectral information can be obtained from Doppler or speckle contrast data. In both the computer simulations and laboratory measurements, we obtain the same spectrum by both methods.

Equation 1, the basis of our analysis, requires the assumption that an ensemble of speckle intensities collected from points distributed in space will have the same variance as an ensemble collected in a time series. This assumption is reasonable for skin, where there is no underlying fixed matrix to modulate the speckle pattern, and hence no fixed speckle pattern, but may not apply in other situations—for example, imaging through the thinned skull of rats in cerebral vessel imaging.^17- 18

Forms of Eq. 1: Display Formula

The typical use of Eq. 1 in speckle contrast analysis has been to choose an appropriate autocorrelation function $Ct(τ)$ with some characteristic time parameter $τc$, based on the assumption that a presumed velocity distribution for the red blood cells (RBCs) will generate an autocorrelation function directly linked to its velocities.^5,19- 20 Substituting the chosen autocorrelation function and then integrating Eq. 1 generates a function relating $τc$ and contrast. Measurement of a single speckle contrast value then allows the characteristic time $τc$ to be measured, and from there an estimate of an average mean RBC velocity can be made. This analysis can be extended to incorporate multiple exposures: Parthasarathy et al. ¹⁸ use multiple-exposure speckle imaging and a speckle model revised from Briers’s to include static speckle patterns, giving an improved linearity in relative flow measurements using flow tubes. Similarly, Smausz et al. ²¹ use multiple exposures and a Lorentzian model for velocity distribution, and add extra fitting parameters to their $K(τc)$ function to obtain a better estimate of $τc$ than those obtained by a single-exposure measurement.

These methods rely on choosing the correct physical model for light scattering in tissue and moving blood, including an appropriate velocity distribution.¹⁹ The point has been made that the statistics of angular scatter from RBCs means that the Doppler spectrum is strongly biased away from simple speed averaging toward lower frequencies because backscatter has a low probability.²² The choice of velocity distribution is also unclear.¹⁹

Once the autocorrelation function and its related power spectrum are measured using multiple exposures and the analysis shown here, speckle measurements are in the same position regarding interpretation as Doppler measurements, and follow the same path to an estimate of blood flow. This interpretation requires a number of assumptions—these are, however, common to all of the analyses here. The assumptions are that (1) there is a fixed matrix of tissue scatterers, with all moving scatterers being RBCs; (2) the velocity distribution of the RBCs is independent of their spatial position; (3) the individual blood cells move independently; and (4) there is an RBC concentration low enough that each scattered photon encounters no more than one RBC.^22- 23 These assumptions will clearly not be met in all cases, particularly when imaging blood vessels directly, but in the case of dermal perfusion, they should be reasonable. The difference between vessel and tissue imaging should also be considered when designing flow tubes in scattering phantom experiments to test perfusion measurements, as used by several groups.^24- 26 Although this approach may provide good calibration for measurement of flow in vessels, the extension of these calibrations to microvascular perfusion measurement may not be reliable.

Bonner and Nossal,²² Nilsson,²³ and other workers in laser Doppler flowmetry give the equation, for an arbitrary velocity distribution,Display Formula

Some problems remain with a direct interpretation of such a perfusion index, whether measured by multiexposure speckle contrast or by Doppler. These problems relate to the assumptions required earlier, particularly the assumption that only the movement in the tissue is blood flow. Some authors report a biological zero, a perfusion index remaining during occlusion of blood flow in Doppler measurements,²⁷ and we have found the same phenomenon in speckle measurements, suggesting that the tissue matrix should not, in fact, be considered totally static. Further measurements of the spectra of occluded tissues, or other biological specimens where there is no net flow, might elucidate the limits of this assumption.

It may be possible, given sufficient knowledge of the expected tissue spectrum, to fit a parametric function to the data at some stages of the laser speckle analysis and then find a perfusion index as a function of the fitting parameters. We have used this approach in skin measurements with some success in previous work,¹² generating a perfusion index that changes with dermal vasodilation and is consistent, on the same subject, over the course of several weeks. This parametric-fitting approach gives up some robustness, particularly the capability to handle any velocity distribution of the moving scatterers, in favor of a simpler measurement requiring fewer exposures.

Multiexposure speckle contrast measurements, using the analysis described, can recover the same spectral information as laser Doppler measurements. Speckle measurements can therefore measure dermal perfusion, with the same assumptions as required for laser Doppler perfusion measurement, but with the advantage over Doppler measurements of generating full images at video rates. This technique can therefore provide an improvement for any current application of laser Doppler imaging.

Multiexposure speckle need not rely on choosing a particular intensity autocorrelation function to generate an estimate of blood flow but can instead calculate such a function from measurements, by analysis of the speckle contrast versus exposure curve.

This work was funded by the Foundation for Research, Science and Technology, Contract No. C08X0201.