1 Introduction

Laser Doppler perfusion imaging (LDPI) [1, 4, 6, 7, 20, 28, 30, 33] and related techniques such as Laser Speckle Contrast Analysis (LASCA) [6–8, 11, 14, 16, 17, 33] and Laser Speckle Imaging (LSI) [9, 26] are established techniques for determining skin perfusion maps, for instance, to diagnose burns [5, 25, 33, 34] or diabetic ulcers [22], to study cerebral blood flow in humans [27] and small animals [13, 15, 24], and for drug uptake studies (e.g. using iontophoresis) [10].

In LDPI, the skin is illuminated with coherent laser light. A fraction of the laser light that entered the tissue interacts with moving red blood cells and obtains a Doppler shift. When these photons leave the tissue they interfere on the detector with each other and with light which has not obtained a Doppler shift, resulting in a time-fluctuating signal. The moments of the power spectrum S(ω) of this photoelectric current are given by

$$ M_i \propto \int\limits_{0}^\infty{\omega^{i} S(\omega)\hbox{d}\omega} $$
(1)

In LDPI, the zeroth order moment (i = 0) is a measure for the concentration of red blood cells whereas the first-order moment (i = 1) is a measure for the flux or perfusion.

In commercially available LDPI devices the area under investigation is scanned with a narrow laser beam. Signal processing of the photoelectric current is performed by analog frequency-weighted filters. Using beam scanning, a perfusion image of 256 × 256 takes approximately 5 min. This long scanning time impedes the observation of fast perfusion changes, for instance, during reperfusion after occlusion. In order to obtain real-time reperfusion images a refresh rate of approximately 25 Hz is needed. Also, for the patient, a short imaging time is desirable (e.g. for burn patients, and in general for young children and elder patients). The need for reduced imaging time is solved by the introduction of CMOS imaging arrays in the field of LDPI [12, 29–31]. This ascent in optical array sensor technology allowed the real-time acquisition of dynamic laser speckle patterns generated by the tissue perfusion, which in principle allows real-time LDPI. Recently, it is shown that the relation between the first moment and tissue blood perfusion still holds for the full-field illumination scheme [12].

In order to take advantage of these developments, the conversion of measured photo current fluctuations into perfusion estimates should be as efficient as possible. A straightforward way to calculate the perfusion in LDPI with a CMOS-camera is to make use of the Fast Fourier Transform (FFT) to calculate the power spectral density (PSD) followed by a weighted integration of the PSD over all the frequencies. Calculating the perfusion using the FFT is time consuming because it has to be taken for each pixel and the spectrum has to be squared and integrated over all the frequencies to obtain one number for the perfusion. Obtaining a perfusion image of 256 × 256 pixels takes up to 10 s on a personal computer with a 3.2 GHz Pentium 4 processor and 1.0 Gb RAM-memory.

Here, we present and evaluate a calculation algorithm for the first-order moment working in the time domain (referred to as TD) which involves less computational steps and requires less data to be transported and stored. The term ‘time domain’ refers to the fact that the perfusion parameter is determined from the time series of photo current fluctuations, without the need to perform a mathematical transformation of this time series, such as a Fourier transform.

2 Method

2.1 Theory

Koelink et al. [23] described a time domain algorithm that estimates the first-order moment of the photo detector current power spectrum via the second-order moment. For this purpose, the assumption was made that the power spectrum is exponentially shaped. In our approach, we directly estimate the first-order moment of the power spectrum.

The convolution property states that the Fourier transform of the product of two functions is the convolution of the Fourier transform of these functions:

$$ \int\limits_{-\infty}^{\infty}{g(t)h(t)\exp(-i\beta t)\hbox{d}t} ={\frac{1}{2\pi}}\int\limits_{-\infty}^{\infty}{G(\omega )H(\beta -\omega)\hbox{d}\omega } $$
(2)

where G(ω) and H(ω) are the Fourier transforms of g(t) and h(t), respectively. For β =  0, this becomes

$$ \int\limits_{-\infty}^\infty{g(t)h(t)\hbox{d}t}={\frac{1}{2\pi}} \int\limits_{-\infty}^\infty{G(\omega)H(-\omega)\hbox{d}\omega} $$
(3)

Taking the time derivative in the time domain, this corresponds with multiplication with iω in the frequency domain:

$$ {\frac{\partial^{n}}{\partial t^{n}}}f(t)\to(i\omega)^n F(\omega ) $$
(4)

Hence, introducing \(h(t)={\frac{\partial}{\partial t}}g^{*}(t)\) and applying the property that the Fourier transform of g *(t) equals G *(− ω), Eq. 3 reduces to

$$ \left\langle g(t){\frac{\partial}{\partial t}}g^{*}(t)\right\rangle ={\frac{1}{2\pi}}\int\limits_{-\infty}^\infty{i\omega G(\omega ) G^*(\omega)\hbox{d}\omega} $$
(5)

We introduce the real signal f(t) (i.e. f *(t) = f(t)) and amplitude spectrum A(ω), which are a Fourier pair:

$$ f(t)={\frac{1}{2\pi}}\int \limits _{-\infty}^{\infty}A(\omega)\exp (i\omega t)\hbox{d}\omega $$
(6)

as well as the power spectrum S(ω) of f(t) which is defined as

$$ S(\omega)\equiv\left|{A(\omega)}\right|^{2}=A^{*} (\omega)A (\omega ) $$
(7)

Equation 5 can then be written as

$$ \left\langle f(t){\frac{\partial}{\partial t}}f (t)\right\rangle = {\frac{1}{2\pi}}\int\limits_{-\infty}^{\infty}{i\omega S(\omega) \hbox{d}\omega } $$
(8)

Because iωS(ω) is an anti-symmetric function, taking the integral from − ∞ to ∞ will result in the right-hand side of Eq. 8 equal to 0. It can also be proven that for a statistically stationary function f(t), the left-hand side of Eq. 8 equals zero [18]. Hence, although correct, Eq. 8 is not useful for our purposes. However, in the field of LDPI the first moment is originally defined as [15]:

$$ M_1\equiv\int\limits_{-\infty}^{\infty}{\left| \omega \right|S(\omega )\hbox{d}\omega } $$
(9)

In Eq. 8, the right-hand side being zero is the result of a balance between the integrals from − ∞ to zero and zero to ∞. Replacing iω with |ω| on the right-hand side of Eq. 8 makes the right hand side unequal to zero, and equal to M 1 as defined by Eq. 9. Equally, we can break the balance in the left-hand side of Eq. 8 by replacing it with the time average of \(\left| f \frac{\partial f}{{\partial t}} \right|.\) Next, we assume that the equal but opposite parts which cancel both the right and left sides of Eq. 8, are equal to each other. This is still a mainly intuitive step, for which mathematical justification is being searched. However, it is a tempting step leading to

$$ M_1\approx 2\pi\left\langle{\left|{f(t){\frac{\partial f(t)}{\partial t}}}\right|} \right\rangle $$
(10)

We emphasize that Eq. 10 will probably provide an approximation of the first-order moment as defined by Eq. 9 rather than being an exact result of straightforward mathematics. Hence, for calculating the first order moment of the power spectrum, we can reduce the number of computational steps in comparison with the FFT-based method : Eq. 10 shows that the first-order moment can be estimated from the fluctuating part of the photocurrent signal and its time derivative. Since

$$ \left| {f(t){\frac{\partial f(t)}{\partial t}}} \right| = \hbox{sign} \left( {f(t){\frac{\partial f(t)}{\partial t}}} \right)f(t){\frac{\partial f(t)}{\partial t}} $$

mathematical justification of Eq. 10 may be found by studying the spectral properties of sign(f(t)∂f(t)/∂t). In our future research, we will investigate the mathematical justification of Eq. 10 along this line of argument.

2.2 Materials

When the light intensity fluctuations are digitally recorded, the input signal f(t) is converted into a discrete signal consisting of samples f i . Equation 10 can be discretized as

$$ M_1={\frac{2\pi\nu_{\rm s}}{N}}\sum\limits_{i = 0}^{N-1}{\left|{f_i \left( {f_{i + 1}-f_i}\right)}\right|} $$
(11)

where νs is the sample frequency, f i the ith sample of f, and N the total number of samples of f. Trials were performed with various discretizations of Eq. 11. Taking the time derivative with f i-1 and f i or f i and f i+1 gave comparable results.

Equation 10 and its discretization in Eq. 11 represent the TD-algorithm and give an expression for the first-order moment of S(ω) in terms of the samples of f. Hence, using the TD-algorithm gives the opportunity to determine the perfusion with some simple mathematical operations in the time domain which do not require large computational power or time.

In order to test the accuracy of the TD-algorithm, several scattering objects were imaged with a homebuilt Laser Doppler Perfusion Imager (LDPI) using a 10-bit CMOS-camera and a 300-mW laser emitting at 671 nm to illuminate an area of approximately 10 × 10 cm. This setup is described in detail in reference [11]. A total of 1,024 images of 128 × 128 pixels were obtained at a frame rate of 27 kHz, resulting in 16,384 raw time traces of 1,024 samples each. Afterwards, DC-subtraction (i.e. subtracting the mean value) was applied to focus on fluctuations in the comparison rather than on the mean value. Without this subtraction, the outcome of the TD-algorithm would mainly be determined by the offset introduced by the DC-value rather than by the fluctuations.

The first sample consisted of a cylindrical piece of Delrin of 40 mm in diameter with a hole of 4 mm in diameter filled with IntraLipid 20% (Fresenius Kabi) and placed on green surgery paper. The second comparison was performed on human tissue. The ring and little fingers of the right hand of a volunteer (male, 26 yr), positioned on a black cloth, were imaged. The last comparison was performed during an occlusion and the reperfusion afterwards [19]. The occlusion was applied by inflating a blood pressure cuff around the upper arm of a healthy subject (male, 27 yr). The measurement on the wrist was started 3 min after the occlusion was applied. Approximately, 5 s after the start of the measurement the occlusion was released. The data was afterwards processed with both algorithms. For every perfusion image of 128 × 32 pixels, the average value was determined based on the complete image, and plotted as a function of time.

Based on these time traces, the first-order spectral moment was calculated with the TD-algorithm and FFT-algorithm. For the FFT-algorithm, we integrated over all the frequencies present in the signal, i.e. from 0 to 13.5 kHz. Calculations were performed using Labview 7.1 (National Instruments) and Matlab 7.2 (the MathWorks) on a personal computer with a 3.2 GHz Pentium 4 processor and 1.0 Gb RAM-memory. For comparison, a scatter plot is made in which the value obtained with the FFT-algorithm is plotted as a function of the value obtained with the TD-algorithm, as well as a linear fit through these values (solid line) and its 95% confidence interval. As guidance to the eye, the line FFT = TD (dotted line) is drawn and the correlation coefficient (r) is given. This allows to analyse the correlation between the two estimates.

Furthermore, the two algorithms are compared using the Bland–Altman method [2, 3]. This involves plotting the difference of the values obtained with the two methods versus their average value. The black solid line in the Bland–Altman plot indicates the mean difference, and the dashed lines indicate the mean plus and minus twice the standard deviation of the difference between the estimates. When the points are normally distributed, 95% of all points lays between these two dashed lines. In order to visualize the trend in the data, the data is divided into horizontal bins with a size of 1,000. In these bins the mean value is determined and plotted (red line).

Each point in the scatterplot as well as in the Bland–Altman plot corresponds with the analysis of the time trace of one pixel.

In order to compare both algorithms with respect to calculation time, the time needed to calculate the perfusion for an image of 256 × 256 pixels for different lengths of the time signal was measured for both algorithms using MatLAB and LabVIEW.

3 Results

In order to give an impression of the variation in spectral content included in our measurements, three types of power spectra are shown in Fig. 1. The spectra have been obtained in regions of 5 × 5 pixels, situated in the image of the Delrin part of the first object, the Intralipid part of the first object, and the tip of the finger.

Fig. 1
figure 1

The power spectra of Delrin (blue line), Intralipid (red line) and human tissue (black line) averaged over an area of 5 × 5 pixels

Full size image

In the case of the cylindrical Delrin sample containing some Intralipid, the resulting images of the first-order moments with FFT and TD algorithms, the scatter plot of FFT-results versus TD-results, and Bland–Altman plot are shown in Fig. 2a–d, respectively. The slope of the linear fit of the scatter plot is 0.8229 with a 95% confidence interval of [0.8221, 0.8238].

Fig. 2
figure 2

128 × 128 pixel images of the perfusion obtained from a cylindrical piece of Delrin of 40 mm in diameter with a hole of 4 mm in diameter, and placed on green surgery paper; the hole was filled with IntraLipid 20% (imaged area approx. 7 × 7 cm). Perfusion plots with a the TD-algorithm, b the FFT-algorithm and c scatter plot of FFT-results versus TD-results (solid line: linear fit; dashed line: TD=FFT) and d Bland–Altman plot of both perfusion plots, the red line indicates the average value of TD-FFT over horizontal intervals of 1,000

Full size image

The results obtained from the hand of the volunteer are shown in Fig. 3. The slope of the linear fit in the scatter plot is 1.1319 with a 95% confidence interval of [1.1314, 1.1324].

Fig. 3
figure 3

128 × 128 pixel images of the perfusion of the finger tips of the hand of a volunteer (imaged area approx. 4 × 4 cm). a Perfusion plot with the TD-algorithm, b perfusion plot with the FFT-algorithm, c the scatter plot of FFT-results versus TD-results (solid line: linear fit; dashed line: TD=FFT) and d Bland–Altman plot of both perfusion plots; the red line indicates the average value of TD-FFT over horizontal intervals of 1,000

Full size image

Figure 4 shows the perfusion during an occlusion and reperfusion calculated with the FFT-algorithm (black line) and TD-algorithm (red line).

Fig. 4
figure 4

Perfusion as a function of time during occlusion and reperfusion calculated with the FFT-algorithm (black line), and the TD-algorithm (red line)

Full size image

The calculation time with both algorithms using MatLAB and LabVIEW as function of the number of raw speckle images involved in the calculation of the spectral first moment, is shown in Fig. 5 along with the calculation time as reported by Serov et al. using the FFT-approach [31].

Fig. 5
figure 5

Calculation time for the TD-algorithm and the FFT-algorithm for an image of 256 × 256 pixels for two software packages, as a function of the number of raw images

Full size image

4 Discussion

4.1 Comparison of both algorithms

The spectra in Fig. 1 show that static Delrin causes the least broadening, while Intralipid shows the highest spectral broadening, as must be expected. Figure 1 gives an impression of the large variety of spectral shapes in the data sets which are used to compare the TD- and FFT-algorithms.

As shown in Figs. 2, 3, the perfusion maps for the several measured samples differ slightly between both algorithms. The scatter plots for the perfusion (Figs. 2c, 3c) show a narrow cloud instead of a straight line indicating that there is some difference between the TD-algorithm and the FFT-algorithm. This difference can also be seen in the Bland–Altman plots (Figs. 2d, 3d) where there is a deviation from the black solid line (i.e. the mean difference), and the cloud is not equally distributed around the solid line. The differences in the perfusion images seem largest for the images of skin perfusion (Fig. 3). However, this is mainly a visual effect due to the large colour differences within the range of image values. A more quantitative comparison is given by the scatter plots and Bland–Altman plots. The agreement between TD and FFT in terms of absence of random variations, is the best for the first-order moment values obtained in static regions, and hence, on cloth material and Delrin. In these regions, first-order moments are found in the interval between 0 and 2.5, where the scatter plots show relatively small random variations compared to the first-order moments exceeding values of 2.5. In these regions, the perfusion image is largely caused by noise in the photo current. In this case, the signal is more evenly spread in the spectral domain than signals caused by Brownian motion or real perfusion. This suggests that the TD-signal can deal better with high- than low-frequency signal components. This difference can have several reasons: a possible source of error in the TD-results is the discretization applied in this algorithm.

In spite of the differences between the TD- and the FFT-algorithm as observed in the scatter plots, the perfusion images obtained with both methods are similar. This similarity can be shown in two ways. First, when a ROI is selected on the finger in Fig. 3a, b (not shown), the standard deviation within these two ROIs is between the 10 and 15%. This is in the same order as the difference between the values of the two ROIs. Therefore, the difference between the TD- and the FFT-algorithms can largely be explained based on statistical grounds. Secondly, the average curve in the BA plots (red line) shows that the systematic difference between TD and FFT is approximately 15%. Note that the large excursions of the mean difference for high image values are in those regions where only a few values contribute to the BA-plots. This is also in the same order as the standard deviation within the two ROIs. Therefore, both methods give similar results and will lead to the same judgment by the physician. Clinical judgment will be based not on separate pixel values, but on the perfusion in a larger tissue area. Therefore, the difference between the FFT- and the TD-algorithms can be seen as noise superimposed on the image plus a systematic difference of maximal 15%. The interpretation by the physician will not be influenced by the difference between the perfusion maps.

The small difference between the TD- and the FFT-algorithms is also visible in the comparison of both algorithms during the occlusion–reperfusion procedure (Fig. 4). With both algorithms, the heart beat is nicely visible after release of the occlusion. The average time between the heart beat peaks is 920 ms, which corresponds with the heart beat frequency of 65 beats per minute. In the reperfusion stage, an offset of < 5% is visible. For low perfusion levels (i.e. during the occlusion), this offset has increased.

In contrast to LDPI, the lower integration limit for the FFT-algorithm was set to 0 Hz, which implies that all the frequencies are taken into account, and both algorithms should be equally sensitive to low frequencies. Noise correction was not applied in both algorithms. This explains why in both algorithms perfusion is still visible during the occlusion.

4.2 Calculation speed

As shown in Fig. 5 the TD-algorithm is at least twice as fast as the FFT-algorithm, so the TD-algorithm is a good alternative for the FFT-algorithm in LDPI. The more the raw images are taken into account, the bigger the saving of time is. Tests (out of the scope of this article, and so not shown), to determine the accuracy of the calculated perfusion level as a function of the number of raw images, reveal that at least 512 images, but preferably more, should be taken into account to obtain a proper perfusion image.

4.3 Implementing TD on high speed CMOS-cameras

CMOS-cameras which are equipped with a programmable memory for data acquisition can implement the discretized version of the TD-algorithm easily which brings (semi-)real-time perfusion imaging within reach. Memory requirements are low, since the TD-algorithm only needs two raw frames for obtaining a preliminary perfusion plot. Every new raw frame can then be used for producing a new perfusion plot by combining it in a weighted way with the previous perfusion plot and the previous raw frame. In this way, the perfusion map will converge to a statistically accurate value, and a perfusion value update can be performed when a new frame is available. In this way, perfusion variations in time can be followed with minimum transport and data storage. The memory requirements are slightly increased by the fact that, in practice, the raw signal must be high pass filtered to suppress the influence of signal artefacts which often occur at frequencies below 20 Hz. One way of implementing high-pass filtering in time domain is by making use of spectral inversion [32]. In this technique, the low-pass filtered signal obtained by a moving average filter is subtracted from the unfiltered signal. Since the low-frequency components are subtracted from the original signal, only the high-frequency components appear in the output, and the signal is high pass filtered. When spectral inversion is implemented in a moving way, it can be incorporated in the TD-algorithm. The use of a moving average approach implies that the number of frames to be stored in memory is larger than two, but still much lower than the amount needed to perform a complete Fast Fourier Transform. In this way, low frequencies can be filtered out, although it will probably not be as good as high-pass filtering in the FFT-algorithm. This means that the time-domain algorithm is possibly more sensitive to low frequencies in the power spectrum, caused, for instance, by heart-beat related rhythmic blood volume variations (the photoplethysmographic signal) [21] or motion artefacts.

5 Conclusion

In conclusion, we have shown that the TD-algorithm can be a good alternative for the FFT-algorithm. The TD-algorithm provides a suitable approximation of the first-order moment rather than being an estimator based on a rigorous mathematical analysis. The new algorithm can determine the perfusion in LDPI twice as fast as the FFT-algorithm and gives equivalent perfusion images. Furthermore, the TD-algorithm gives the possibility to perform on-the-fly estimation because it only needs a few new frames to update a previous perfusion estimation. This results in lower memory requirements compared to the memory needed to perform an FFT, which is a major advantage of the TD-algorithm.