Dynamic speckle imaging of human skin vasculature with a high-speed camera

Simon Erdmann; Simon Erdmann; Flora Weissgerber; Élise Colin Koeniguer; Xavier Orlik

doi:10.1364/OE.452111

1. Introduction

The idea of using the statistics of a dynamic speckle pattern to infer parameters of the movement of light scatterers is at least 40 years old [1–3]. Most of the work since then has been based on extracting parameters of the motion from the autocorrelation function of the signal. In particular, the so-called correlation time of the speckle signal appears crucial. It corresponds both to the time lag needed for the autocorrelation to fall below a certain value and to the characteristic time taken by the scatterers to change their layout significantly. Yet, its main flaw resides in the difficulty to provide a universal definition of this latter notion: it depends on features of the imaging system (numerical aperture, illumination wavelength) but also on the regime of the motion of the scatterers. Some authors suggest extracting a characteristic distance from the system properties and then considering the ratio of this distance to the correlation time as the speed of the scatterers. Yet, there is no consensus about the characteristic length which fits best this purpose [4]. On the other hand, assuming a specific regime for the scatterers’ motion, a rigorous link can sometimes be made between correlation time and a physical parameter. For instance, some models link correlation time with the diffusion coefficient $D$ of a Brownian motion [5], or with the standard deviation $\delta v$ of the speed distribution of a ballistic motion [5].

However, for decades, such an assumption was then difficult to validate. Recording a moving speckle signal at a sufficient acquisition rate and with an exposure time low enough to sample the autocorrelation function correctly was made impossible by the available hardware, at least in full field. Therefore, motion monitoring from speckle images has widely relied on speckle contrast, defined as the ratio of the standard deviation to the mean of the signal intensity. The more the speckle decorrelates of itself during the exposure time, the more the recorded signal appears as blurred, which mathematically translates into a loss in contrast. Through a relatively simple general integral expression, the contrast can be expressed as a function of the autocorrelation function [5], and, as a result, of the correlation time of the signal. For practical applications, a given shape of the autocorrelation function, corresponding to a scatterers motion model, can be used in this expression. Numerous such models have been suggested, all adapted to a certain type of motion: the Lorentzian model [2] is designed for Brownian motion, the Gaussian model [6] for ballistic motion, the rigid-body model [4] for scatterers moving together, etc. Alternatively, if the exposure time is large enough compared to the correlation time, the ratio between the latter and the former can be estimated with $1/C^2$, where $C$ is the speckle contrast [7]. This estimation is efficient in its validity range [8,9], even though it only provides approximate results.

Yet, regardless of the exact method used, contrast-based speckle imaging has always remained qualitative and quite imprecise. With these methods, autocorrelation models needed to have only one parameter to be invertible, and thus were necessarily overly simple. Moreover, there are numerous cases where the motion of the scatterers is not known well enough to make a physically well-grounded choice of model. Thus, the model chosen was often unreliable and impossible to validate. Some techniques, such as Multi-Exposure Speckle Imaging [10,11] have worked around these issues by repeating the acquisition with several exposure times, but at the price of a longer acquisition process and with still a limited accuracy. On the other hand, autocorrelators or fast single-pixel detectors enabled to access the short-term dynamics of scatterers, through autocorrelation or zero-crossings [3], but only for one point at a time. There seemed to be a gap between single-pixel techniques, based on expensive hardware, like Doppler imaging or the ones mentioned above, and laser speckle imaging, which was affordable and able to provide high-resolution images, but remained qualitative [12]. At the end of the 2000s, the first full-field Doppler imaging enabled real-time vasculature imaging, but at the price of a still coarse spatial resolution [13–15].

However, in recent years, cameras with high temporal resolution have become widespread, allowing the recording of data with an acquisition frequency high enough to sample the autocorrelation function itself correctly. Fitting the empirical autocorrelation to a given model is consequently becoming a leading technique to estimate movement [16–18]. This enables both to refine models by adding new parameters and to validate them by inspecting the quality of the fit.

Our work is part of this general trend. With the help of a high-speed camera (up to several dozen kHz of acquisition rate), we tried to directly image the skin vasculature of a human subject. Yet, we chose not to rely on a particular physical model. Our purpose was to explore the abilities of a simple high-speed imaging setup when used on various body parts. This involved being agnostic concerning the shape of our autocorrelation functions.

Whereas the abilities of speckle vasculature imaging were thoroughly studied for animal models [11,18] or the human retina [8,19], the intact human skin was comparatively little studied. This can be explained by the fact that the skin vasculature is not as easy to see as the brain or retina vessels. Moreover, keeping a human subject still might be more difficult than an anesthetized mouse.

Yet, there are many applications where qualitative imaging of the skin would be a considerable help to physicians in seeing if vasculature is present. For instance, monitoring its development is essential in cancerology (as abnormal growth may be the sign of a tumor) or in skin grafting (as monitoring the recovery of blood vessels is essential after a skin graft).

Consequently, in this paper, our purpose is to image the skin vasculature of a human subject, keeping both the hardware setup and the algorithms simple. Concerning the hardware, the subject will not be constrained in any way, but will only do his best to stay still. As for the method in itself, we will compare simple imaging parameters (autocorrelation coefficient, number of zero-crossings), adapted to a high-speed setup, with the speckle contrast, better fit for long integration times, for reference. Our purpose is to lay the foundations of skin vasculature imaging methods making the most of a high-speed camera. This means identifying what imaging parameter fits best this purpose, how to process the data to make the best of it, and the constraints and limits of such a setup in terms of acquisition parameters. In particular, making the most out of the high-speed camera means trying to minimize total acquisition time. We point out that this time is also the amount of time during which a potential patient has to stay still: reducing it may increase their comfort and limit the need for post-processing to correct unintentional motion. Moreover, this may enable to access non-stationary phenomena such as the heartbeat. In the longer term, we also want the proposed method to be ready for real-time imaging.

2. Materials and methods

2.1 Available data

We used a similar protocol to image different body parts of the same 50-year-old Caucasian male subject.

For each acquisition, the skin was illuminated from a typical distance of $\sim$50 centimeters with a near-infrared LASER (RGB Photonics Lambda Mini @ 785 nm). The linearly polarized LASER beam was made slightly divergent by removing its internal convergent lens. It enabled us to reach a compromise between spatial power density and size of the illuminated area. The output power of the LASER was 100 mW, and the spatial power density over the samples was about 30 W/m².

We point out that we chose not to enforce the spatial Nyquist-Shannon criterion, which requires that speckles are at least twice as large as pixels. Whereas it would have increased the contrast of the raw signal, it would also have implied a smaller aperture, which, in turn, would have reduced the illumination power received by the sample, to a point our camera was not sensitive enough for. Note that, even though many teams take care of the Nyquist-Shannon criterion, the theory does not strictly require it. Our work is then is then a chance to prove experimentally that enforcing it is not necessary.

Then, the illuminated area was filmed with a Phantom VEO 710L camera from about $\sim$50 centimeters, with an acquisition frequency of 2000 to 30,000 Hz and spatial sampling of 128$\times$128 to 320$\times$240. The number of raw images acquired was constrained by the internal memory of the camera, such that the maximum size of an acquisition file was about $\sim$500 MB. This corresponds to a total of about 50 million pixels (10 bits/pixel), and ultimately to a few seconds of total acquisition time.

Concerning the choice of the body parts imaged, our purpose was to test our setup on various vasculature scales and biological phenomena. The six datasets we present here try to reflect this diversity:

• The first three datasets were acquired by imaging the anterior surface of the left wrist, with respective acquisition frequencies of 3000, 2000, and 10,000 Hz. These datasets will be named WRIST_1, WRIST_2, and WRIST_3. As the wrist contains major veins and arteries, we expect to detect rather large vessels.
• The fourth and fifth datasets respectively consist of images of the end of two fingers of the left hand (dataset FINGERS) and of the right ear (EAR). They were respectively acquired at 30,000 and 2000 Hz. Both fingers and ears are vascularized by small vessels.
• The last dataset corresponds to the palm of the left hand, which was gently scraped beforehand to trigger a superficial inflammation in the shape of a smiley face. It was then filmed with an acquisition frequency of 10,000 Hz. We do not expect to retrieve vasculature from this dataset but to image the inflammation, as it is invisible on the raw images of the stack.

For each acquisition, the subject remained sitting on a chair. For the acquisitions involving the hand or the wrist, the hand was simply lying on the optical table, with the palm pointing to the camera. For the dataset EAR, the head was not maintained and the subject only did his best to stay still.

The main acquisition parameters are summarized in Table 1.

Table 1. Main acquisition parameters of the datasets used for this study.

View Table

All these data constitute Dataset 1 [20], which is publicly available. An example image of each stack is shown in Fig. 1.

Fig. 1. Example raw image for each acquisition

Download Full Size | PDF

2.2 Preprocessing

2.2.1 Image registration

Despite the very short acquisition times, it proved impossible for the filmed subject to remain completely still compared to the system resolution during the acquisitions.

Therefore, the images were not perfectly coregistered, as we could note by manually examining different frames successively. This turned out not to be an issue for FINGERS, WRIST_3, and PALM, for which the displacement during the acquisition was small enough. For the other three datasets (WRIST_1, WRIST_2 and EAR), this was a problem. Consequently, we decided to investigate how to coregister the images.

2.2.1.1 Coregistration algorithm

Regarding the coregistration algorithm itself, we chose not to reimplement it ourselves and to use available code, namely the registration algorithm offered by OpenCV, through its Python interface [21].

This algorithm is based on a gradient descent on a similarity criterion, like the Lucas-Kanade algorithm [22]. The main advantage of gradient-descent methods is that they provide by construction a sub-pixel shift value, whereas correlation-based algorithms with an exhaustive search are discrete. Yet, as they are local, gradient-descent algorithms are unable to retrieve large shifts and tend to fall into local minima. This issue is solved by performing hierarchical motion estimation, that is by estimating the shift at different scales successively. This feature is also implemented in OpenCV using a Gaussian pyramid.

2.2.1.2 High-pass spatial filtering for wrist datasets

The raw images of the EAR dataset have enough visible features to be coregistered using a conventional algorithm. In the wrist datasets, however, the only structures we can hope to lean on are a few hairs in the upper right corner of the image.

Yet, if we directly apply our coregistration algorithm to these images, it fails to retrieve the shifts, according to a simple visual inspection. This is likely due to the hairs not being prominent enough, as the algorithm seems to estimate the flow based on a more visible structure: the laser spot. This is due to the fact that we could not make the sample illumination completely uniform: this would have required a much more divergent beam. This, in turn, would have provided too little spatial power density, and caused quantization noise and shot noise. As a result, the laser spot is visible on raw images. The problem is that it remains stationary, as the laser does not move with the subject. Therefore, based on the immobility of the spot, the algorithm estimates a zero flow on the whole image.

To solve this difficulty, we noted that the laser spot appears as a slow-varying gradient on the images, whereas the hair is thin structures, consequently corresponding to high spatial frequencies. That is why we removed the low-frequency components of the images, using a Gaussian high-pass filter of size 4 $\times$ 4, to suppress the non-uniformity of the illumination while preserving the hair structures. Note that this spatial filtering was only performed to compute the shifts between images, and was not applied during preprocessing.

Eventually, the coregistration process can be summarized as follows: for every raw image, 1. a spatial Gaussian high-pass filter of size 4 $\times$ 4 was applied, 2. the image was coregistered with the first image of the stack, thanks to the algorithm provided by OpenCV. We point out that our coregistration process relies on visible features and is not specially designed for speckle images. This process is adapted to a specific situation (the need to cut the low frequencies due to the laser spot), and required the specific features we found in the images (ear edges and wrist hair). Automating the process would require supplementary work, as well as adapting it to situations where no prominent features are visible. However, we think the details we gave may still be used by others. More generally, we believe coregistration could often be useful in dynamic speckle methods.

Visual inspection let us think our coregistration process was successful. The benefits of coregistration are more thoroughly developed in Section 3.5.

2.2.2 High-pass temporal filtering

Most of the processing we perform on the images theoretically requires that the temporal signals are stationary. On the other hand, despite coregistration, the subject’s movement tends to generate sudden local changes in the mean of the temporal signals. We believe this is because a given point of the imaged object is illuminated differently over time, as the illumination is not completely uniform and is independent of the subject’s motion. This, in turn, creates artefacts on the final images.

To keep the signal average constant for each pixel, we applied a temporal Gaussian high-pass filter with a standard deviation corresponding to 32 ms (for instance, 64 images at a sampling rate of 2000 Hz). The Gaussian high-pass filter $\mathcal {G}^{hp}$ is defined by:

\mathcal{G}^{hp} = \mathcal{I} - \mathcal{G}^{lp}

where $\mathcal {I}$ is the identity and $\mathcal {G}^{lp}$ is the standard Gaussian low-pass filter. Thus, the Gaussian high-pass is complementary to the Gaussian low-pass: the former cuts the frequencies the latter keeps and vice versa.

The usefulness of the whole preprocessing scheme (coregistration and temporal filtering) is evaluated at the end of the Results Section, in Sec. 3.5.

2.3 Methods to image vascular activity

2.3.1 Temporal contrast

Temporal contrast has been used for a long time to map activity from speckle images. It works well if the integration time is significantly larger than the decorrelation time of the speckle pattern, i.e. if the speckle becomes independent to itself several times during integration.

However, in our configuration, the integration time is very short. Therefore, we can expect it to be shorter than the decorrelation time. That is why we calculated the temporal contrast not only on raw images but also after simulating a longer integration time: assuming that the sensor is perfectly linear, we can simulate an exposure time $N$ times longer by adding pixel values temporally along blocks of $N$ contiguous images.

In all cases, the contrast is computed on the temporal axis in order to maximize spatial resolution. Considering the set of values taken by a pixel over time as our statistical sample, empirical contrast was estimated by dividing the sample standard deviation by the sample mean.

2.3.2 First term of the autocorrelation function

In many recent works, motion parameters are extracted from the empirical autocorrelation function [17,18,23]. However, in most of these works, the authors choose to fit a model to their empirical autocorrelation function to extract parameters. These models can sometimes be difficult to justify physically, such as the single or double exponent model [17,23], or quite complex and only suitable to a particular physical situation [18]. Moreover, our autocorrelation curves appear to be too noisy after a few terms to allow for a proper fit. In addition, our goal remains to obtain a result in quasi real-time, which prohibits any computationally expensive processing.

Therefore, we chose a simpler method and directly used the first term of the autocorrelation function, i.e. the temporal correlation coefficient between consecutive images, as our base imaging parameter. The $n-$th term of the autocorrelation function is defined by:

(1)$$r_n = \frac{{\mathbb{E}\!\left[I(t + n\Delta t) I(t)\right]} - {\mathbb{E}\!\left[I(t\right]}^2}{{\mathbb{V}\!\left[I(t)\right]}}$$

where:

• $\Delta t$ is the acquisition sampling step, that is the reciprocal of the acquisition sampling rate;
• As the signal is stationary, $t = m\Delta t$ with $m \in \mathbb {N}$ is any positive multiple of $\Delta t$ such that $t + n\Delta t < N \Delta t$, where $N$ is the length of the discrete signal;
• $I(t)$ is the temporal signal (i.e. the time-integrated intensity) of the pixel of interest at discrete time $t$;
• $\mathbb {E}$ and $\mathbb {V}$ respectively refer to the expected value and the variance, which we effectively estimate with their canonical estimators.

$r_1$ corresponds to the case where $n=1$.

Our imaging parameter of choice will be $1 - r_1$: as $r_0 = 1$, $1 - r_1$ is proportional to the opposite of $(r_1 - r_0) / \Delta t$, which is the discretization of the derivative of $r$ at time lag $0$. In other words, $1 - r_1$ is a short-term decorrelation rate. As a result, we consider it a measurement of the activity, or motion, within the considered pixel.

We point out that many authors rather use $g_2$ instead of $r$ as their autocorrelation function. The discrete version of $g_2$ is defined by:

(2)$$g_2[n] = \frac{{\mathbb{E}\!\left[I(t + n \Delta t) I(t)\right]} }{{\mathbb{E}\!\left[I(t)\right]}^2}$$

Therefore, $r_n$ and $g_2[n]$ are related through the equations:

(3)$$\begin{aligned} r_n & = \frac{{\mathbb{E}\!\left[I(t)\right]}^2}{\mathbb{V}\!\left[{I(t)}\right]} (g_2[n] - 1)\\ & = \frac{g_2[n] - 1}{C^2} \end{aligned}$$

where $C$ is the speckle contrast. Inverting this equation provides:

$$g_2[n] = 1 + C^2 r_n$$

So the three quantities are linked through a simple relation.

$r$ typically decreases from 1 to 0, whereas $g_2[n]$ decreases from $1 + C^2$ to 1. Therefore, as $C$ is not constant across the image, the range of $g_2[n]$ is not normalized, contrary to $r_n$. $g_2[n]$ is consequently less comparable between two pixels.

Another perspective on this point results from rewriting (3):

\begin{aligned} r_n &= \frac{g_2[n] - 1}{C^2} \\ &= \frac{g_2[n] - 1}{g_2[0] - 1} \\ &= \frac{g_2[n] - g_2[+\infty]}{g_2[0] - g_2[+\infty]} \\ \end{aligned}

This emphasizes the idea that $r_n$ can be considered a normalized version of $g_2[n]$. We considered it made $r_1$ more adapted to our method than $g_2[1]$. For example, a $r_1$ value of $0.9$ means the same (that is a strong short-term correlation) for any pixel, whereas a $g_2[1]$ value of $1.1$ cannot be interpreted without also knowing the contrast at the pixel considered (one cannot know whether it is close to $1 + C^2$, and consequently if the correlation is high or not).

2.3.3 Link between contrast and autocorrelation

2.3.3.1 Siegert relation

Before developing the links between autocorrelation and contrast, we shall mention the importance of the Siegert relation in dynamic speckle imaging. The Siegert relation links the continuous autocorrelation function $g_2$ of the intensity over a pixel to the continuous autocorrelation function $g_1$ of the field complex amplitude. $g_1$ and $g_2$ are respectively defined by:

\left\{ {\begin{array}{l} {{g_1}(\tau ) = \frac{{{\mathbb{E}}\left[ {{{\mathbf A}^*}(t){\mathbf A}(t + \tau )} \right]}}{{{\mathbb{E}}\left[ {|{\mathbf A}(t){|^2}} \right]}}\quad \quad (4)}\\ {{g_2}(\tau ) = \frac{{{\mathbb{E}}\left[ {{I_c}(t){I_c}(t + \tau )} \right]}}{{{\mathbb{E}}{{\left[ {I(t)} \right]}^2}}}n\quad \quad (5)} \end{array}} \right.

where $I_c(t)$ is the continuous time intensity, such that $I(t) = \int _{t-T/2}^{t+T/2} I_c(u) \,\mathrm {d}u$.

With this in mind, the Siegert relation states that:

(6)$$g_2(\tau) = 1 + \beta {|g_1(\tau)|}^2$$

where $\beta$ is a parameter comprised between $0$ and $1$. $\beta$ measures the contribution to the total signal coherence of factors other than time, such as space ($\beta$ decreases when there are several speckles per pixel) or polarization (for example, $\beta$ increases if the incident wave is linearly polarized, or if one linear polarization is selected out of the backscattered wave). $\beta$ is also the squared contrast of the underlying continuous intensity signal over one pixel:

(7)$$\beta = \frac{ {\mathbb{V}\!\left[I_c(t)\right]} }{ {{\mathbb{E}\!\left[I_c(t)\right]}}^2 }$$

A demonstration of this result is provided as Supplement 1 (Sec. S1). This section of Supplement 1 also includes a further discussion about $\beta$. Its detailed expression in a particular case (when only loss of spatial coherence is taken into account).

Note that a consequence of the Siegert relation is that the continuous autocorrelation coefficient $r$ of $I_c$ is equal to ${|g_1|}^2$:

\begin{aligned} r(\tau) &= \frac{{\mathbb{E}\!\left[I_c(t) I_c(t+\tau)\right]} - {{\mathbb{E}\!\left[I_c(t)\right]}}^2}{ {\mathbb{V}\!\left[I_c(t)\right]} } \\ &= \frac{{{\mathbb{E}\!\left[I_c(t)\right]}}^2}{{\mathbb{V}\!\left[I_c(t)\right]}} \frac{{\mathbb{E}\!\left[I_c(t) I_c(t+\tau)\right]}}{{{\mathbb{E}\!\left[I_c(t)\right]}}^2} - \frac{{{\mathbb{E}\!\left[I_c(t)\right]}}^2}{{\mathbb{V}\!\left[I_c(t)\right]}} \\ &= \frac{1}{\beta} g_2(\tau) - \frac{1}{\beta} \end{aligned}

where we used the definition (5) of $g_2$ and (7) for $\beta$. Inverting (6) then brings $r(\tau ) = {|g_1(\tau )|}^2$.

2.3.3.2 Contrast and autocorrelation

the more the speckle pattern decorrelates of itself during the exposure time, the more the contrast is reduced. Therefore, $C$ is a priori expected to be an increasing function of the correlation time $\tau$ when the exposure time is fixed. This principle is formalized by the integral formula linking $C$ to $r$, the continuous autocorrelation function [5]:

(8)$$C^2 = \frac{2\beta}{T} \int_{0}^T \left[ 1 - \frac{\tau}{T} \right] r(\tau) \,\mathrm{d} \tau$$

where $T$ is the exposure time, and $\beta$ is the parameter appearing in the Siegert relation. The contrast actually measures the loss of coherence of the discretized signal compared to the continuous field. With this in mind, the integral term in (8) corresponds to the loss of time coherence due to time integration, whereas $\beta$ stands for the other sources of coherence loss, as stated previously.

For given $\beta$ and $T$, we therefore expect $C$ to increase with $\tau$. As we will assume that the continuous autocorrelation $r$ is decreasing from 1 to 0, it can easily be shown that $r_1$ should, as well, be an increasing function of $\tau$. Therefore, we expect $C$ to be high where $r_1$ is high and low where $r_1$ is low.

Yet, this might depend on the precise statistical behavior of the signal at the examined point, as we will show by comparing the spectral behaviors of both parameters.

2.3.3.3 Spectral behavior

For this section, we will examine one single pixel. Let $s$ be its time signal, and $S$ the Discrete Fourier Transform (DFT) of $s$. $S$ is defined by:

(9)$$S[k] = \sum_{n=0}^{N-1} e^{{-}i\frac{2kn\pi}{N}} s[n]$$

With these notations, the squared contrast $C^2$ at the considered pixel is given by:

\begin{aligned} C^2 &= \frac{\frac{1}{N} \sum_{n=0}^{N-1} s^2[n] - {\left( \frac{1}{N} \sum_{n=0}^{N-1} s[n] \right)}^2 }{{\left( \frac{1}{N} \sum_{n=0}^{N-1} s[n] \right)}^2} \\ &= \frac{\frac{1}{N} \sum_{n=0}^{N-1} s^2[n]}{{\left( \frac{1}{N} \sum_{n=0}^{N-1} s[n] \right)}^2} - 1 \end{aligned}

The sum in the denominator is $S[0]$, by definition of $S$. As for the numerator, the Parseval equality provides:

\sum_{n=0}^{N-1} |s[n]|^2 = \frac{1}{N} \sum_{k=0}^{N-1} |S[k]|^2

This equality is easily demonstrated by injecting (9) into the right-hand side and then developing it. Then, we have:

\begin{aligned} C^2 &= \frac{\frac{1}{N^2} \sum_{k=0}^{N-1} {|S[k]|}^2}{\frac{1}{N^2} {S[0]}^2} - 1\\ &= \frac{\sum_{k=0}^{N-1} {|S[k]|}^2}{{S[0]}^2} - 1 \end{aligned}

So $1/C^2$ is given by:

\frac{1}{C^2} = \frac{1}{\frac{\sum_{k=0}^{N-1} {|S[k]|}^2}{{S[0]}^2} - 1}

${|S[k]|}^2$ is the power spectrum of the signal, and the sum in the numerator is its total energy. This means that $1/C^2$ is an increasing function of $\frac {{S[0]}^2}{\sum _{k=0}^{N-1} {|S[k]|}^2}$, that is the ratio of the signal energy contained in its mean, i.e. in its DC component.

On the other hand, we will show $1 - r_1$ can be interpreted as the spectral spread, or width, of this energy.

First, (1) provides:

1 - r[1] = \frac{R[0] - R[1]}{R[0] - {\left( \frac{1}{N} \sum_{n=0}^{N-1} s[n] \right)}^2}

where $R$ is the non-normalized autocorrelation function of the signal:

(10)$$R[\Delta n] = \frac{1}{N} \sum_{n=0}^{N-1} s^*[n] s[n+\Delta n]$$

where the signal was assumed to be periodic (i.e. $s[N + n] = s[n]$). This assumption is not physically grounded, but will have little consequence as we will only use it for $n=1$ in a sum of $N$ terms. It will make the calculations simpler and clearer.

Now, the Discrete Fourier Transform of $R$ is $1/N$ times the power spectrum of the signal, that is $\frac {1}{N}|S|^2$. This result can be shown by computing the Discrete Fourier Transform of (10), or considered resulting from the definition of the power spectrum of a discrete random signal [24]. It is also analogous to the Wiener-Khinchin Theorem for deterministic functions [25]. Therefore, we can write $R[0]$ and $R[1]$ as inverse Fourier transforms:

\begin{aligned} 1 - r[1] &= \frac{ \frac{1}{N^2} \sum_{k=0}^{N-1} {|S[k]|}^2 - \frac{1}{N^2} \sum_{k=0}^{N-1} e^{i\frac{2k\pi}{N}} {|S[k]|}^2 }{ \frac{1}{N^2} \sum_{k=0}^{N-1} {|S[k]|}^2 - \frac{1}{N^2} {|S[0]|}^2 } \\ &= \frac{ \sum_{k=0}^{N-1} \left( 1 - e^{i\frac{2k\pi}{N}} \right) {|S[k]|}^2 }{ \sum_{k=1}^{N-1} {|S[k]|}^2 } \\ \end{aligned}

As ${|S[k]|}^2$ is the discrete Fourier transform of a real signal (namely $R$), it is symmetrical with respect to the Nyquist frequency, which corresponds to $k=N/2$. Therefore, $S[k] = S[N-k]$ for every $k$. Consequently:

\begin{aligned} \sum_{k=0}^{N-1} \left( 1 - e^{i\frac{2k\pi}{N}} \right) {|S[k]|}^2 &= \sum_{k=1}^{\left\lfloor \frac{N-1}{2} \right\rfloor} \left( 1 - e^{i\frac{2k\pi}{N}} \right) {|S[k]|}^2 + 2 \epsilon_N {|S[N/2]|}^2 \\ &\qquad + \sum_{k=1}^{\left\lfloor \frac{N-1}{2} \right\rfloor} \left( 1 - e^{i\frac{2(N-k)\pi}{N}} \right) {|S[k]|}^2 \end{aligned}

where $\epsilon _N = 1$ if $N$ is even and $0$ if $N$ is odd, and the term for $k=0$ is not shown because it is zero.

\begin{aligned} \sum_{k=0}^{N-1} \left( 1 - e^{i\frac{2k\pi}{N}} \right) {|S[k]|}^2 &= \sum_{k=0}^{\left\lfloor \frac{N-1}{2} \right\rfloor} \left[ (1 - e^{i\frac{2k\pi}{N}}) + (1 - e^{i\frac{2(N - k)\pi}{N}}) \right] {|S[k]|}^2 + 2 \epsilon_N {|S[N/2]|}^2 \\ &= \sum_{k=0}^{\left\lfloor \frac{N-1}{2} \right\rfloor} \left[ 2 - e^{i\frac{2k\pi}{N}} - e^{{-}i\frac{2k\pi}{N}} \right] {|S[k]|}^2 + 2 \epsilon_N {|S[N/2]|}^2 \\ &= \sum_{k=0}^{\left\lfloor \frac{N-1}{2} \right\rfloor} \left[ 2 - 2\cos\left( \frac{2k\pi}{N} \right) \right] {|S[k]|}^2 + 2 \epsilon_N {|S[N/2]|}^2 \\ &= \sum_{k=0}^{\left\lfloor \frac{N-1}{2} \right\rfloor} 4\sin^2\left( \frac{k\pi}{N} \right) {|S[k]|}^2 + 2 \epsilon_N {|S[N/2]|}^2 \\ &= 2 \sum _{k={-}\left\lfloor \frac{N - 1}{2} \right\rfloor} ^{\left\lfloor \frac{N}{2} \right\rfloor} \sin^2\left( \frac{k\pi}{N} \right) {|S[k]|}^2 \\ \end{aligned}

$1 - r[1]$ is then:

(11)$$1 - r[1] = \frac{ 2 \sum_{k={-}\lfloor (N-1)/2 \rfloor}^{\lfloor N/2 \rfloor} \sin^2 \left( \frac{k\pi}{N} \right) {|S[k]|}^2 }{ \sum_{k=1}^{N-1} {|S[k]|}^2 }$$

The $\sin$ function is increasing on the interval $[0, \pi /2]$. Therefore, $k \mapsto \sin ^2 \left ( \frac {k\pi }{N} \right )$ is an increasing function of $k$ on $\{ 0, \ldots, \lfloor N/2 \rfloor \}$. So according to (11), $1 - r_1$ can actually be interpreted as the spread, or the width, of the power spectral density of the signal. This makes a $1 - r_1$-based method close to full-field Doppler imaging [13–15], but at a lower computational price (as it does not require to compute the DFT of the signal).

Therefore, from a spectral point of view, $1 - r_1$ and $1/C^2$ measure two different features of the power spectrum: the former is sensitive to the spectral width of the signal, that is to its high-frequency components, whereas the latter is sensitive to its DC component. Therefore, the fluctuations of the underlying continuous signal will tend to contribute to $1/C^2$ if they are noticeably faster than the sampling rate, and to $1 - r_1$ if they are slower.

2.3.4 Zero-crossings

The zero-crossings method consists of counting the number of times a signal of interest crosses its mean (and consequently, if it is centered, the number of times it crosses zero) in a given time interval. If the signal is discrete, as is our speckle signal for one pixel, what is counted is the number of times the sign of the centered signal changes between two consecutive points.

The method was applied on several occasions to laser speckle signals [3,26]. It relies on thorough theoretical studies of the average number of zero-crossings for various kinds of processes [27–29]. In particular, the autocorrelation function of a process has been linked to its average number of zero-crossings in various cases. If $X$ is a discrete stationary Gaussian signal of mean 0 and length $N$, it may be shown that:

(12)$$\frac{{\mathbb{E}\!\left[N_0\right]}}{N-1} = \frac{\arccos r_1}{\pi}$$

where $N_0$ is the random variable counting the zero-crossings of $X$. A demonstration of this formula is proposed in Supplement 1 (Sec. S2). As $\arccos$ is a decreasing function, it means that ${\mathbb {E}\!\left [N_0\right ]}$ can be interpreted as a spectral width of the signal as well as $1 - r_1$. Moreover, if $\tilde {N}_0$ is an estimator of $N_0$ (typically, the count of the actual number of sign changes), this formula can be inverted to provide an alternative estimator of $r_1$:

(13)$$\tilde{r}_1' = \cos\left(\frac{\pi\tilde{N}_0}{N-1} \right)$$

3. Results

3.1 Comparison between activity retrieved from temporal contrast and from autocorrelation

As we explained sooner, since the acquisition time is expected to be shorter than the correlation time of the speckle pattern, speckle contrast on raw stacks is likely to give poor results. Therefore, for comparison between contrast and autocorrelation, we computed:

• The first term of the autocorrelation function $r_1$;
• The first term of the autocorrelation function $g_2$;
• The temporal contrast on raw stacks;
• The temporal contrast with a synthetic exposure time of 1 ms, obtained by assuming the camera response is linear and averaging consecutive images by blocks of an appropriate size;
• The temporal contrast with a synthetic exposure time of 10 ms, obtained with the same method as for 1 ms.

The “raw stacks” used for each of these operations were preprocessed beforehand, as explained in Section 2.2.

$r_1$ was then used to obtain the imaging parameter $1 - r_1$, and we computed $1/C^2$ from the contrast. As for $g_2[1]$, we chose to directly image it as it was. The resulting images are shown in Fig. 2, together with the temporal mean of each stack for reference. Note that the colormap is inverted for $g_2[1]$, to keep the same color order as in the other images.

Fig. 2. Different imaging parameters computed for each dataset. The six columns respectively show the temporal mean (after coregistration) of each stack, $1-r_1$, $g_2[1]$, $1/C^2$ on raw stacks, $1/C^2$ after synthesizing an exposure time of 1 ms, and $1/C^2$ after synthesizing an exposure time of 10 ms. On the $1 - r_1$ image of PALM (the second image of the bottom row), we added red arrows which point to the eyes and nose of the smiley face.

Download Full Size | PDF

Our observations are the following:

• On the $1 - r_1$ images, presented in the first column, we can see different linear structures, that we interpret as blood vessels of various sizes. On the PALM dataset, $1 - r_1$ reveals the smiley face drawn by scraping the skin. Its nose and eyes are pointed by red arrows on the image. The only dataset for which no structure is clearly visible is WRIST_3. We can also notice that some parts of these structures can be guessed on the temporal mean images of some datasets, especially those of the wrist. This means that the imaged vessels have a slightly different absorption coefficient from the skin, which is consistent with our expectation of large vessels in the wrist. Yet, the considerable gain in visibility on the $1 - r_1$ images suggests the information it carries does not only come from absorption.
• Compared to $1 - r_1$, the $g_2[1]$ images might show similar structures on some of the datasets, but with a very low contrast. Remember that $g_2[1]$ is a mixture of $r_1$ and $C$, according to the formula: $$$g_2[1] = 1 + C^2 r_1$$$ Thus, taking $r_1$ removes the contrast information from $g_2[1]$ by normalizing it. This normalization appears strongly beneficial. However, to our knowledge, $g_2[1]$ is never used as an imaging parameter anyway. It only seemed important to us to have it appear because $g_2$ is often used in the literature as part of a model where the effect of the contrast is accounted for.
• As for contrast images, the $1/C^2$ images at 1 ms exhibit roughly similar structures as $1 - r_1$ on most stacks. Yet, similarly to $g_2$, their resolution is visually much coarser and fewer details are visible.
• The $1/C^2$ images at 10 ms are very similar to $1/C^2$ images at 1 ms, but look noisier. This is easily explained by the fact the contrast was estimated on fewer images (as the time windows corresponding to one image are longer).
• On the $1/C^2$ images at full frequency, results are even worse: vascular structures are invisible for WRIST_1, WRIST_3 and FINGERS, and barely visible for WRIST_2 and EAR.

We also point out that our choice to not enforce the Nyquist-Shannon criterion did not prevent us from obtaining high-resolution $1 - r_1$ images. This is an important point, as this criterion used to be considered a mandatory prerequisite to apply historical speckle imaging methods [30].

3.1.1 Differences between $r_1$ and the contrast

As emphasized above, $1/C^2$ gives better results after simulating an exposure time of 1 ms. This confirms that the contrast is an efficient parameter only if the exposure time is long enough.

We also noticed differences in quality between $1/C^2$ and $1 - r_1$ images. Apart from them, several discrepancies between these images are unexpected. As explained in the Methods section, $C$ should roughly be an increasing function of $r_1$ — and therefore $1/C^2$ should also be an increasing function of $1 - r_1$. Yet, this principle is not always respected. For instance, let us look back at our images obtained from the EAR dataset, reproduced in Fig. 3(a) and 3(b). There are some areas (annotated on these same figures) where evolutions of $1 - r_1$ and $1/C^2$ are contradictory, that is where $1 - r_1$ is high and $1/C^2$ is low or the opposite.

Fig. 3. Various parameters for dataset EAR. Areas of interest, where $1/C^2$ and $1 - r_1$ do not vary together as expected, are annotated on subfigures (a) and (b).

Download Full Size | PDF

These areas are:

1 The concha (central hole of the ear): high $1-r_1$, low $1/C^2$;
2 The scapha (top left of the auricle): high $1-r_1$, low $1/C^2$;
3 The apparent convergence area of vessels: high $1-r_1$, low $1/C^2$;
4 The fossa triangularis (top right of the auricle): low $1-r_1$ around the vessel, high $1/C^2$.

A general explanation of these differences may be provided by (8). This equation expresses the fact that $C$ decreases when temporal coherence decreases (as formalized by the integral term including $r$) but also with all other factors likely to reduce coherence for other reasons (loss of spatial coherence, polarization, etc.), which are all contained in $\beta$. These factors have no reason to have a constant influence throughout the image. Therefore, a possible explanation for the signal behavior in the regions of interest is that the loss in coherence is relatively low in regions 1, 2, 3 (where $1/C^2$ is low, so $C$ is high) and high in region 4. For example, more volume scattering may lead to a higher loss of coherence through depolarization. As for spatial coherence, it can be influenced by the roughness of the filmed object, through the spatial frequencies of its bumps and holes. All these phenomena may influence $\beta$.

We also have in mind two other explanations, which are justified by the particular features of the signal. Let us take a look at the mean and standard deviation of the dataset. Figure 3(c) and 3(d) show the temporal mean and standard deviation of the signal for every pixel. The temporal speckle contrast is therefore a ratio of these two images. We notice that our regions of interest are either particularly dark (1, 2, 3) or bright (4) in terms of the signal mean, and consequently of light intensity.

First, the darker areas (1, 2, 3) are possibly much noisier than the rest of the image: due to shot noise, the signal-to-noise ratio (SNR) of a camera is commonly an increasing function of the camera response. Under this hypothesis, some white noise would have been added to the signal in darker regions, increasing contrast but also reducing autocorrelation, as such noise is self-uncorrelated.

Second, another explanation could justify the behavior of region 4. Throughout our study, we have always assumed the camera response was perfectly linear, that is the value of a pixel is proportional to the light intensity it received during the exposure time. Yet, it is known that most cameras are not naturally linear: they typically are for the lower response values but tend to saturate at higher values [31]. Therefore, the local derivative of the camera response as a function of intensity tends to decrease. This involves that, at high response values, the computed speckle contrast can be artificially reduced compared to lower values. This phenomenon might explain the behavior of region 4, which could be verified by properly linearizing the camera.

3.2 Comparison with zero-crossings

As detailed in Section 2.3.4, (12) provides a direct formula expressing the zero-crossings count as a function of $r_1$. Therefore, we do not expect it to give us any new piece of information compared to the $1 - r_1$ images. In order to validate the formula, we computed a second estimator of $r_1$ from zero-crossings based on (13) for acquisition WRIST_1. We call this estimator $r_1^{(zc)}$, while the canonical estimator of $r_1$ is called $r_1^{(direct)}$. Then, we plotted the points $(r_1^{(direct)}, r_1^{(zc)})$ and computed the parameters of the linear regression (without intercept) between the estimators. The plot can be seen in Fig. 4(a), and the parameters of the linear regression were $a = 1.003$ (slope) and $R^2 = 0.948$ (square of the linear correlation coefficient). Therefore, for the kind of precision that we need, both estimators are equivalent. This fact is visually confirmed by Fig. 4(b) and 4(c), where $r_1$ and $N_0$ images are placed next to each other: as $N_0$ is a monotonous function of $r_1$, both images exhibit roughly the same features.

Fig. 4. Comparison between direct $r_1$ and $N_0$ for dataset WRIST_1.

Download Full Size | PDF

To make an informed choice concerning our estimator of $r_1$, we performed Monte-Carlo simulations (not detailed here) which tended to show the direct estimator has a slightly (but probably not decisively) smaller variance. Yet, the zero-crossings parameter may still be interesting in some specific cases, such as situations where computational efficiency is critical and zero-crossings turn out to be significantly faster to compute.

3.3 Influence of acquisition frequency and exposure time

As mentioned previously, WRIST_3 is an exception among our datasets, since its $1 - r_1$ image only exhibits barely visible structures. On the other hand, some structures are visible on images of our other datasets of the wrist (WRIST_1 and WRIST_2). These datasets also have a lower acquisition frequency: respectively 3000 Hz for WRIST_1 and 2000 Hz for WRIST_2, whereas it is 10,000 Hz for WRIST_3.

Therefore, we decided to simulate acquisition parameters closer to those of the first two wrist acquisitions, that is a longer exposure time, by averaging the stack by blocks of consecutive images. In order to more generally study the effect of the acquisition rate and the exposure time on the $1 - r_1$ images, we did the same with every other stack. We obtained the results shown in Fig. 5.

Fig. 5. $1 - r_1$ with various acquisition frequencies. Each column corresponds to a given frequency, except for column 3, for which $\Delta t = 1/2000$ or $1/3000$ second depending on datasets. This is due to the synthetic acquisition frequency having to be a divisor of the base frequency, what constrained its possible values.

Download Full Size | PDF

We can see that, with an acquisition frequency of 2000 Hz or less, the $1 - r_1$ image of WRIST_3 lets appear vascular structures. On PALM, taking $\Delta t = 1 / 2000$ s also increases the visibility of the smiley face. On the other hand, it is interesting to notice that, with other acquisitions, decreasing the acquisition frequency degrades the image quality. Even with WRIST_3 and PALM, the image quality decreases below 2000 Hz.

The first result means that, with our method, there is no benefit to increasing the acquisition frequency beyond a few kHz to image the wrist, as we cannot take advantage of the gain in time resolution. It is interesting to note that this is not the case with the fingers: despite an acquisition frequency of 30 kHz, finger vessels were visible on the $1 - r_1$ image of FINGERS (Fig. 2(d), left column). According to (11), which provides a spectral interpretation of $r_1$, this means that the signal has much more high-frequency content in FINGERS than in WRIST_3. Intuitively, this should mean that scatterers tend to move faster in the fingers than in the wrist.

Thus we saw that a too high acquisition frequency may degrade the imaging quality of $1 - r_1$. In addition, we notice that using a low acquisition frequency also gives poor results. Therefore, we conclude that there is an optimum range for the acquisition frequency for $1 - r_1$, out of which the imaging quality is poor. This range depends on the body part of interest: in particular, it is much higher for fingers than for wrists.

3.4 Minimum total time needed to obtain an activity image

In previous sections, for each dataset, we used the whole stack to compute imaging parameters. As one of our goals is to generate a final image as fast as possible in terms of total acquisition time, we tried to determine how low we could have it while keeping the vascular structures distinguishable from the background.

For each dataset, we extracted subsets of consecutive images from the full raw stack, corresponding to respective time intervals of 500, 50, and 5 milliseconds. Then, we coregistered the images, applied our temporal high-pass filter, and computed $1-r_1$ (after simulating an acquisition frequency of 2000 Hz for WRIST_3). Results are shown in Fig. 6.

Fig. 6. $1 - r_1$ on substacks for each datasets. Columns respectively correspond to time intervals of 500, 50 and 5 ms.

Download Full Size | PDF

For all datasets, the minimum time needed to obtain a readable image seems to be around 50 ms — maybe a bit less for WRIST_1 and a bit more for the others. This means a good-quality activity image can be obtained in at most a few hundreds of milliseconds. This is very encouraging both in the perspective of making the acquisition process instantaneous for a patient or a physician, and of displaying the real-time evolution of vasculature at a good refresh rate (around 10 Hz).

3.5 Benefits of preprocessing

As mentioned in the Methods section, all our datasets were preprocessed, first with coregistration for some of them, then with high-pass temporal filtering. As this processing is not standard in the field of dynamic speckle, we will demonstrate its relevance here.

For this purpose, we computed $1-r_1$ with or without coregistration, with or without filtering, on datasets WRIST_2 and EAR. Results are shown in Fig. 7.

Fig. 7. $1 - r_1$ with different preprocessing schemes on datasets WRIST_2 and EAR.

Download Full Size | PDF

We can make the following observations:

• Without coregistration nor filtering, wrist hair appears on the $1-r_1$ image of WRIST_2. They disappear after coregistration, but vertical artifacts tend to appear. These artifacts vanish after filtering. When images are filtered, coregistration also seems to improve visual resolution.
• With the EAR dataset, no vascular structure at all appears without filtering. When filtering, coregistration also helps increase resolution, like with WRIST_2.

These results are the most spectacular ones among our datasets. Yet, we have seen no case where coregistering or filtering degrades the quality of $1-r_1$ images. Therefore, it may be a good idea to get used to systematically perform these two operations. If not, one should keep in mind that on datasets where $r_1$ looks unable to recover any information, like EAR, high-pass filtering can make a big difference.

4. Conclusion

In this study, we first tried to compare the results of the dynamic speckle technique obtained on human skin with three different parameters: the contrast, the first term of autocorrelation, and the zero-crossings count.

The most obvious noteworthy result we got is that the zero-crossings method and the first term of the autocorrelation function give the same information. This could be expected due to the formula provided by (12), but, to our knowledge, this formula had never been experimentally verified. We can now state that both parameters provide two good estimators for the same quantity. Using the first term of the autocorrelation directly is more natural, as it is closer to usual statistical processing. However, zero-crossings might be worthy to consider in computationally critical applications.

We also compared the first term of the autocorrelation function, $r_1$, and the contrast $C$. First, a theoretical study enabled us to state that both quantities can measure activity but through different features of the power spectrum of the signal. More precisely, $1 - r_1$ and $1/C^2$ can respectively be interpreted as its spread and the relative weight of its DC component. If the signal is sampled at a sufficient rate, what is the case with high-speed cameras like ours, its fluctuations should tend to contribute to $1 - r_1$. On the other hand, if the acquisition rate is slow and the exposure time is long compared to the correlation time, fluctuations are time-integrated and contribute to the DC component of the signal, and thus to $1/C^2$. Therefore, we could expect that $1 - r_1$ would give better results than $1/C^2$. This statement was confirmed experimentally: whereas we could globally see the same structures in $1 - r_1$ as in $1/C^2$, these structures were systematically more visible in $1 - r_1$ images — both in their extent and in their distinguishability from the background. Moreover, as contrast requires a longer exposure time, data meant to compute $r_1$ are intrinsically faster to acquire than for $C$. We also showed that the acquisition time could be made even shorter: for some datasets, a few dozen raw images are sufficient for this purpose, which correspond to a few dozen milliseconds of acquisition. All this makes $r_1$ a parameter of choice for qualitative motion mapping. Furthermore, among all possible normalizations of the autocorrelation function, our results tend to show $r$ is a proper choice: using $g_2[1]$, for example, provided poor results.

However, our visual results would not have looked the same without our preprocessing scheme. Image coregistration, as expected, and high-pass filtering, maybe more surprisingly, turn out to improve the quality of our $1-r_1$ images substantially. It would be interesting to try to apply the same operations with other dynamic speckle methods, for example implying model fitting on the autocorrelation function.

Ultimately, $r_1$ appears as a simple, robust, and effective imaging parameter for vasculature or skin inflammation imaging. It requires no model or particular assumptions, works with a limited amount of data, and is fast and easy to compute.

Funding

Institut National Du Cancer; Office National d'études et de Recherches Aérospatiales; Région Occitanie Pyrénées-Méditerranée.

Acknowledgments

Simon Erdmann’s thesis work is supported by a grant co-financed by the Occitanie Region and Onera. This work was also financially supported by the ESCAPADS contract, in the framework of the 2014–2019 Cancer Plan rolled out by the French National Cancer Institute (INCa).

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are publicly available in Dataset 1, Ref. [20].

Supplemental document

See Supplement 1 for supporting content.

References

1. A. F. Fercher, “Velocity measurement by first order statistics of time-differentiated laser speckles,” Opt. Commun. 33(2), 129–135 (1980). [CrossRef]

2. A. F. Fercher and J. D. Briers, “Flow visualization by means of single-exposure speckle photography,” Opt. Commun. 37(5), 326–330 (1981). [CrossRef]

3. N. Takai, T. Iwai, and T. Asakura, “Real-time velocity measurement for a diffuse object using zero-crossings of laser speckle,” J. Opt. Soc. Am. 70(4), 450–455 (1980). [CrossRef]

4. D. D. Duncan, S. J. Kirkpatrick, and J. C. Gladish, “What is the proper statistical model for laser speckle flowmetry?” Proc. SPIE 6855, 685502 (2008). [CrossRef]

5. R. Bandyopadhyay, A. S. Gittings, S. S. Suh, P. K. Dixon, and D. J. Durian, “Speckle-visibility spectroscopy: A tool to study time-varying dynamics,” Rev. Sci. Instrum. 76(9), 093110 (2005). [CrossRef]

6. J. D. Briers and A. F. Fercher,Laser Speckle Technique For The Visualization Of Retinal Blood Flow, in Max Born Centenary Conf, vol. 0369 (International Society for Optics and Photonics, 1983), pp. 22–28.

7. J. W. Goodman, Statistical Optics (John Wiley & Sons, 2000), chap. 6.1, pp. 228–232, 2nd ed.

8. H. Cheng and T. Q. Duong, “Simplified laser-speckle-imaging analysis method and its application to retinal blood flow imaging,” Opt. Lett. 32(15), 2188–2190 (2007). [CrossRef]

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

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

11. A. B. Parthasarathy, S. M. S. Kazmi, and A. K. Dunn, “Quantitative imaging of ischemic stroke through thinned skull in mice with Multi Exposure Speckle Imaging,” Biomed. Opt. Express 1(1), 246–259 (2010). [CrossRef]

12. D. Briers, D. D. Duncan, E. Hirst, S. J. Kirkpatrick, M. Larsson, W. Steenbergen, T. Stromberg, and O. B. Thompson, “Laser speckle contrast imaging: Theoretical and practical limitations,” J. Biomed. Opt. 18(6), 066018 (2013). [CrossRef]

13. A. Serov, B. Steinacher, and T. Lasser, “Full-field laser Doppler perfusion imaging and monitoring with an intelligent CMOS camera,” Opt. Express 13(10), 3681–3689 (2005). [CrossRef]

14. M. Draijer, E. Hondebrink, T. van Leeuwen, and W. Steenbergen, “Twente Optical Perfusion Camera: System overview and performance for video rate laser Doppler perfusion imaging,” Opt. Express 17(5), 3211–3225 (2009). [CrossRef]

15. M. Leutenegger, E. Martin-Williams, P. Harbi, T. Thacher, W. Raffoul, M. André, A. Lopez, P. Lasser, and T. Lasser, “Real-time full field laser Doppler imaging,” Biomed. Opt. Express 2(6), 1470–1477 (2011). [CrossRef]

16. Z. Hajjarian and S. K. Nadkarni, “Evaluating the Viscoelastic Properties of Tissue from Laser Speckle Fluctuations,” Sci. Rep. 2(1), 316 (2012). [CrossRef]

17. M. M. Qureshi, J. Brake, H.-J. Jeon, H. Ruan, Y. Liu, A. M. Safi, T. J. Eom, C. Yang, and E. Chung, “In vivo study of optical speckle decorrelation time across depths in the mouse brain,” Biomed. Opt. Express 8(11), 4855–4864 (2017). [CrossRef]

18. D. D. Postnov, J. Tang, S. E. Erdener, K. Kılıç, and D. A. Boas, “Dynamic light scattering imaging,” Sci. Adv. 6(45), eabc4628 (2020). [CrossRef]

19. Y. Lu and R. K. Wang, “Removing dynamic distortions from laser speckle flowgraphy using Eigen-decomposition and spatial filtering,” J. Biophotonics 15(1), e202100294 (2022). [CrossRef]

20. X. Orlik, S. Erdmann, F. Weissgerber, and E. C. Koeniguer, “Speckle image stacks acquired on human skin with a high-speed camera,” zenodo (2021), https://doi.org/10.5281/zenodo.5802352.

21. G. Bradski, “The OpenCV Library,” Dr. Dobb’s Journal of Software Tools (2000).

22. B. D. Lucas and T. Kanade, “An Iterative Image Registration Technique with an Application to Stereo Vision,” in Proceedings of the 7th International Joint Conference on Artificial Intelligence, vol. 2 (Morgan Kaufmann Publishers Inc., San Francisco, CA, USA, 1981), pp. 674–679.

23. M. Golberg, R. Califa, J. Garcia, and Z. Zalevsky, “Analyzing the requirements of high-speed camera parameters for enhanced laser speckle sensing of flow dynamics,” Engin. Res. Exp. 2(3), 035032 (2020). [CrossRef]

24. A. Papoulis and S. U. Pillai, Probability, random variables and stochastic processes (McGraw-Hill, 2001).

25. D. C. Champeney, A Handbook of Fourier Theorems (Cambridge University Press, 1987), chap. 11, pp. 102–103, 1st ed.

26. H. Cheng, Q. Luo, Z. Wang, L. Yao, S. S. Ulyanov, E. I. Zakharova, and S. Zeng, “Application of laser speckle interferometry for monitoring the dynamics of lymph flow,” in International Workshop on Photonics and Imaging in Biology and Medicine, vol. 4536 (International Society for Optics and Photonics, 2002), pp. 130–136.

27. S. O. Rice, “Mathematical analysis of random noise,” The Bell Syst. Tech. J. 24(1), 46–156 (1945). [CrossRef]

28. N. D. Ylvisaker, “The Expected Number of Zeros of a Stationary Gaussian Process,” Ann. Math. Stat. 36(3), 1043–1046 (1965). [CrossRef]

29. B. Kedem, “On frequency detection by zero-crossings,” Signal Process. 10(3), 303–306 (1986). [CrossRef]

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

31. J. E. Garcia, A. G. Dyer, A. D. Greentree, G. Spring, and P. A. Wilksch, “Linearisation of RGB Camera Responses for Quantitative Image Analysis of Visible and UV Photography: A Comparison of Two Techniques,” PLoS One 8(11), e79534 (2013). [CrossRef]

Name	Description
Dataset 1	Raw image stacks used in the manuscript
Supplement 1	Demonstration of the Siegert relation and the zero-crossings formula

Dataset name	Field of view (cm)	Image size	Sampling rate (Hz)	Exposure time ( $μ$ s)	Number of images	Total duration (s)
WRIST_1	$\sim$ 3	320 $\times$ 240	3000	250	4981	1.66
WRIST_2	$\sim$ 3	256 $\times$ 256	2000	490	5858	2.93
WRIST_3	$\sim$ 3	128 $\times$ 128	10,000	32.9	21,482	2.15
FINGERS	$\sim$ 2	128 $\times$ 128	30,000	32.9	21,482	0.72
EAR	$\sim$ 4	256 $\times$ 256	2000	490	5858	2.93
PALM	$\sim$ 4	256 $\times$ 256	10,000	99,6	5858	0.59

Dynamic speckle imaging of human skin vasculature with a high-speed camera

Abstract

1. Introduction

2. Materials and methods

2.1 Available data

2.2 Preprocessing

2.2.1 Image registration

2.2.1.1 Coregistration algorithm

2.2.1.2 High-pass spatial filtering for wrist datasets

2.2.2 High-pass temporal filtering

2.3 Methods to image vascular activity

2.3.1 Temporal contrast

2.3.2 First term of the autocorrelation function

2.3.3 Link between contrast and autocorrelation

2.3.3.1 Siegert relation

2.3.3.2 Contrast and autocorrelation

2.3.3.3 Spectral behavior

2.3.4 Zero-crossings

3. Results

3.1 Comparison between activity retrieved from temporal contrast and from autocorrelation

3.1.1 Differences between $r_1$ and the contrast

3.2 Comparison with zero-crossings

3.3 Influence of acquisition frequency and exposure time

3.4 Minimum total time needed to obtain an activity image

3.5 Benefits of preprocessing

4. Conclusion

Funding

Acknowledgments

Disclosures

Data availability

Supplemental document

References

Supplementary Material (2)

Data availability

Cited By

Figures (7)

Tables (1)

Equations (25)

Optics Express