Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Correction of overexposure in laser speckle contrast imaging

Open Access Open Access

Abstract

Laser speckle contrast imaging (LSCI) is a method to visualize and quantify tissue perfusion and blood flow. A common flaw in LSCI variants is their sensitivity to the optical setup parameters and that they operate well only on statistics of undistorted laser speckle patterns. The signal saturation of the sensors makes the contrast calculation misleading; hence the illumination level must be well controlled. We describe the theoretical explanation for the saturation-caused degradation. We introduce a linear extrapolation method to eliminate the overexposure induced error up to an extent of 60-70% saturated pixel count. This, depending on the contrast value and use case, enables to use 3-8 times higher external illumination level with no deterioration of the contrast calculation and thus the measured blood flow index. Our method enables a higher signal-to-noise ratio in darker areas by allowing the use of higher illumination, utilizing a larger portion of the dynamic range of the sensors, and making the illumination level setting less cumbersome.

© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

Laser Speckle Contrast Imaging (LSCI) was first introduced for the study of blood perfusion in [1]. The technique is applied in a wide spectrum of areas and used in the biomedical field for non-invasive measurement of the microcirculation and blood perfusion on the brain surface [2] and skin [3], it is used in retina diagnostics [4] and general clinical applications [5,6], as well as for semen motility qualification [7]. LSCI largely builds on the coherent property of laser light; interference occurs when coherent light is reflected from various geometries. An imaging device captures the interference in a from of a random pattern. Then, statistics of the random pattern can be used to derive relative flow speeds. The relation between the laser speckle statistics and flow rate (blood flow index) depends on multiple factors, including the quality of the laser, the optical arrangement, the characteristics of the sensor and the calibration methods [815].

Two of the issues are the over- and underexposure. Underexposure deteriorate the calculated statistics of the speckle pattern and decrease the signal-to-noise ratio due to small discrete values. This issue is addressed in [16] with a correction term. Despite that overexposure alters significantly the speckle statistics as well, for the best of our knowledge, no publication has dealt with the saturation effect. As a rule of thumb, it is recommended to set the average of the detected intensity pattern near or below the middle of the dynamic range of the sensor, so that saturation and clipping is avoided [17]. We argue that even small percentages of the saturated pixels can alter the speckle statistics and derived blood flow indices. Furthermore, the constant need for adjusting the exposure places an additional burden on researchers and developers of LSCI. We believe, that making LSCI methods more robust against saturation, several low intensity and low signal to noise ratio related artifact [18] can be reduced by allowing higher illumination.

In this study, we analyze the effect of overexposure on speckle statistics and propose a compensation method, which provides consistent results in a wide range of exposures. The correction algorithm uses extrapolation of the function of contrast degradation and the ratio of saturated pixels. We show the capabilities of our algorithm with numerical simulations and on real world measurements as well.

2. Laser speckle statistics

LSCI utilizes the statistics of the random speckle pattern based on a measure of contrast $\kappa$. As random pattern changes in accordance with the movement of the scatterers, the flow rate $\propto 1\tau _c$. The different $\tau _c$ decorrelation times are then mapped to a contrast map calculated as:

$$\kappa^2_s = \frac{\sigma^2(I_s)}{\langle I_s \rangle^2},$$
where $\sigma (I_s)^2$ and $\langle I_s \rangle$ are the sample variance and the sample mean calculated in a sliding window of a small neighborhood (e.g $N\times N=5 \times 5$). The reduced second moment $\upsilon _2$ depending on $T$ exposure time and $\tau _c$ is given in the following model [19]:
$$\upsilon_2(T,\tau_c) = \frac{\beta}{T^2} \int_{0}^T \!\!\!\! \int_{0}^T \rho \lvert g_1^n(t'{-}t^{\prime\prime}) \rvert ^2 + (1-\rho) \lvert g_1^n(t'{-}t^{\prime\prime}) \rvert ^2 dt'dt^{\prime\prime},$$
where $\rho$ is the ratio of the dynamic scatterers and $(1-\rho )$ is the ratio of the static scatterers to total scatterers, and $t'{-}t''$ give the time differences during the exposure, $\beta$ is a normalization factor depending on the speckle size [20], detector physical resolution and additional factors of the experimental setup [21,22]. $g_1$ is the electric field autocorrelation, which according to [1] gives us $g_1^n(t) = exp(-(t/\tau _c)^n)$ relation, where $\tau _c$ is the correlation time. The value of $n$ depends on the type of flow dynamics [13]. The contrast values are usually averaged on multiple consecutive frames to improve the signal-to-noise ratio. There are other approaches using temporal or spatiotemporal contrast calculations [10,23,24].

It is important to note, that the outlined methods rely on undistorted speckle patterns, meaning there should be no significant alteration found compared to the ideal Gaussian distribution. Overexposure results in clipped signals and truncated distribution, making the calculations misleading. Our aim is to compensate the saturation effect on contrast calculation during the preprocessing of overexposed speckle patterns.

2.1 Effect of signal saturation

Several previous studies [21,22,25] detail the statistical properties of the speckle pattern. Working with biological samples and realistic sensors, the recorded speckle pattern becomes a composite of multiple origins. Three major factors can be distinguished in this composite: a static element, a dynamic element and a factor, which accounts for the camera noise. The static component can be modeled as a single fully developed speckle pattern. The component of interest is the dynamic one, which is described by the integration of the varying speckle pattern in the time domain over a single exposure. The general composite of the factors can be modeled by the sum of three random variables as [26]:

$$I_s = (1-\rho) S(\lambda_S) + \rho D(\alpha_D,\beta_D) + W,$$
where $S$ is a random variable of exponential distribution, representing static scatterers. Dynamic scatterers are described with the random variable $D$ of gamma distribution. The gamma distribution can be defined by $\alpha _D > 0$ shape and $\beta _D > 0$ rate parameters. Camera noise is taken into account as a random variable of normal distribution $W$ and $\rho$ is the ratio of the dynamic scatterers to total scatterers (as in Eq. (2)).

In the followings, we focus only on large signal behavior, thus we simplify by assuming the camera noise to be negligible $W \approx 0$. Since the sum of different multiple gamma distributions can be estimated by a single near gamma distributed gamma-series [27,28], we handle $I_s$ as a single one, regardless of static scatterer content. The described model does not include several realistic image acquisition and laser source error phenomena (such as laser bandwidth, pixel cross talk, shoot noise, speckle oversampling). In order to handle these effects, we implemented a numerical approach that is described in the following sections.

Let $I_{i}$ be the individual samples with gamma distribution, so the probability density function can be described as:

$$f(I_{i})=\frac{\beta^\alpha}{\Gamma(\alpha)}I_{i}^{\alpha-1}e^{-\beta I_{i}},$$
and the cumulative distribution function is:
$$F(I_{i})=\frac{\gamma(\alpha, \beta I_i)}{\Gamma(\alpha)} ,$$
where $\Gamma (\alpha )=\int _{0}^{\infty }t^{\alpha -1}e^{-t}dt$ is the gamma function, and $\gamma (\alpha,\beta I_i)=\int _{0}^{\beta I_i}t^{\alpha -1}e^{-t}dt$ is the lower incomplete gamma function. In order to calculate the contrast, the expected value and the variance of Eq. (4) is needed. The two values are calculated as $E(f)=\alpha /\beta$ and $\sigma ^2(f)=\alpha /\beta ^2$ based on the contrast calculation equations $\kappa _s^2 = 1/\alpha$ for the unsaturated reference distribution.

Next, let us investigate the saturated case by defining it with the following cumulative (truncated) distribution function:

$$F_{sat}(I_{i})=\begin{cases} F(I_i), &\textrm{if}\; I_i < I_{sat} \\ 1, &\textrm{otherwise}. \end{cases},$$
where $I_{sat}$ is the saturation threshold. The expected value of the truncated gamma distribution in the undistorted $[0, I_{sat}]$ region is:
$$E(\gamma(\alpha, \beta I_{sat})) = \frac{1}{\beta} \frac{\gamma(\alpha+1, \beta I_{sat})}{\gamma(\alpha, \beta I_{sat})}$$

The expected value of the saturated distribution $E(f_{sat})$ can be calculated by weighting the truncated gamma distribution by $F_{sat}(I_{i})$ and adding the saturation threshold $I_{sat}$ weighted by $1-F_{sat}(I_i)$:

$$E(f_{sat}) = \frac{\alpha}{ \beta} \frac{\gamma(\alpha+1, \beta I_{sat})}{\Gamma ( \alpha + 1)} + I_{sat} \left( 1- F_{sat}(I_{i}) \right)$$

Variance is calculated based on the expression of moments $\sigma ^2(x) = E(x^2) - E(x)^2$, and that the combined variance is the sum of the individual variances of independent samples:

$$\sigma^2(f_{sat}) = \frac{\alpha}{\beta^2} \frac{\gamma(\alpha+2, \beta I_{sat})}{\Gamma( \alpha +1)} - E(f_{sat})^2 + I_{sat}^2 \big( 1- F_{sat}( I_{i}) \big).$$

The contrast equals to $1/\alpha$ in the limiting case when $I_{sat} \to \infty$, and it becomes zero when $I_{sat}\to 0$.

For completeness, it is important to mention that the quantization of very low intensity values cause a different kind of distortion in the observed contrast values. As it is described in [16], the effect of the bit depth $b$ on the contrast can be modeled by defining $\kappa ^2$ corrected contrast:

$$\kappa^2 \approx \kappa_s^2 \left( 1 + \frac{I_{sat}}{E(f) (2^b - 1)} \right).$$

Equations (8)–10 enable to simulate and analyze both the effect of under- and overexposure.

2.2 Analytical solution for the reference contrast

We refer to the contrast, which we could observe under ideal conditions as the reference contrast. The reference contrast $1/\alpha$ can be calculated by using higher-order partial derivatives of the saturation ratio and sample mean with respect to the saturation threshold ($I_{sat}$). For ease of notation we introduce saturation ratio $R$ and notation $Y$ as:

\begin{gather}R = 1-F_{sat}(I_{i}) \end{gather}
\begin{gather}Y = E(f_{sat})-I_{sat}R = \frac{\alpha}{ \beta} \frac{\gamma(\alpha+1, \beta I_{sat})}{\Gamma ( \alpha + 1)} \end{gather}

Note, that the saturation ratio is equal to the complementary cumulative distribution function (tail distribution). The shape and rate parameters of the distribution can be calculated in the following way:

\begin{gather}\text{Let} ~~~ A = \frac{\partial^2 Y}{\partial I_{sat}^2} = \frac{1}{\Gamma(\alpha)} \beta^{\alpha} e^{-\beta I_{sat}} \Big[ \alpha I_{sat}^{\alpha-1} - I_{sat}^{\alpha} \beta \Big], \end{gather}
\begin{gather}\text{and let} ~~~ B = \frac{\partial R}{\partial I_{sat}} ={-}\frac{1}{\Gamma(\alpha)} \beta^{\alpha} I_{sat}^{\alpha-1} e^{-\beta I_{sat}}, \end{gather}
\begin{gather}\text{then} ~~~ \frac{A}{B} ={-}\alpha + \beta I_{sat}, \end{gather}
\begin{gather}\frac{\partial}{\partial I_{sat}} \frac{A}{B} = \beta, \end{gather}
\begin{gather}\Big(\frac{\partial}{\partial I_{sat}}\frac{A}{B} \Big) I_{sat}-\frac{A}{B} = \alpha \end{gather}

Though the result is theoretically feasible, it is unsuitable for practical application. First and second order derivatives and their ratio are required for the implementation, which makes the solution sensitive to noise and numerically unstable in the case of small differences.

3. Simulation

In order to investigate the practical effects of saturation, we performed a thorough image-based numerical simulations using the methods and models presented in Refs. [29], [30], [31].

To model a realistic camera and laser source, the following effects are incorporated in the simulations: (i) fixed pattern noise and shot noise modeled by Poisson distribution, (ii) finite bit length data representation by rounding to 12-bit integer, (iii) lateral pixel crosstalk implemented as 7x7 Gaussian kernel, (iv) speckle size, and (v) static scatter content. The simulated wavelength was $\lambda = 820 nm$ and the sensor size was 256x256 pixels. The images were scaled with different illumination values, rounded, then normalized to unity by $I_{sat} = 2^{12}$. The saturation effect was modeled by clipping the intensity.

First, speckle patterns were generated with increasing illumination ($\beta$ rate parameter). The used parameters were $\alpha = 4$, the mean speckle size was 4 pixels, shot noise $\sigma _n = 0.01$, lateral pixel crosstalk $\sigma _p = 0.5$, and $\sigma _\lambda = 3 nm$. First, the illumination level ($I_{unsat}$) was set to the maximum level, where the speckle pattern had no saturated value ($R=0$). Relative to this value, the intensity was varied in $0.25I_{unsat}$, $1I_{unsat}$, $2.5I_{unsat}$, and $4I_{unsat}$ steps. The resulting images and their probability distributions are presented in Fig. 1.

 figure: Fig. 1.

Fig. 1. Density functions of 256x256 pixels of simulated speckle patterns are shown. The used parameters are $\alpha = 4$, speckle size 4 pixels, shot noise $\sigma _n = 0.01$, lateral pixel crosstalk $\sigma _p = 0.5$, static scatterer $\rho =0$. Subplot b) shows the speckle pattern at maximum light intensity ($I_{unsat}$) where no saturation occurred. The subplots a), c), d) correspond to $0.25, 2.5, 4$ times this intensity. The contrast and saturation ratio of the four speckle patterns are $K_s = 0.46, 0.47, 0.41, 0.27$ and $R_s = 0\%, 0\%, 9.1\%, 43.5\%$, respectively.

Download Full Size | PDF

The detailed relation between the contrast, the observed mean illumination and the saturation ratio is shown in Fig. 2. To calculate the saturation ratio $R_s$, the intensity levels are first binarized; 0 is assigned if there is no saturation, 1 if the pixel value reached $I_{sat}$. Then, the binary map is averaged in $N\times N$ neighborhoods:

$$\begin{aligned} R_s = \langle bin(I_i) \rangle \text{, where} \end{aligned}$$
$$bin(I_i) = \begin{cases} 0, & \textrm{if}\; I_i < I_{sat} \\ 1, &\textrm{otherwise}. \end{cases}$$

 figure: Fig. 2.

Fig. 2. Subplot a) shows the contrast values and the saturation ratio as a function of the mean observed illumination. The yellow markers corresponds to the numerical simulation of Fig. 1. Subplots b) and c) show the linear and the exponential components of the contrast approximation as a function of the saturation ratio.

Download Full Size | PDF

It can be seen that the contrast value decreases as more and more pixels become saturated. The reason behind this is that the mean value of the pattern is increasing asymptotically to the saturation level, while the standard deviation is decreasing.

The contrast $\kappa _s$ can be approximated as a function of the saturation ratio $R_s$, with a linear and two exponential components:

$$\kappa_s(R_s) \approx p_1 \cdot R_s + p_0 + a_1e^{{-}b_1R_s} - a_2e^{{-}b_2(1-R_s)}.$$
Fig. 2(c) shows the detrended, exponential component only, with contrast values for the same examples.

Next, we investigated the contrast deterioration as a function of various flow rates and imaging parameters: (i) as a function of $\alpha$ shape parameter, (ii) speckle size, (iii) pixel cross talk and (iv) shot noise. It can be seen in Fig. 3 that the main tendency of contrast deterioration does not alter under various scenarios.

 figure: Fig. 3.

Fig. 3. The upper subplots show the contrast value as a function of the saturation ratio under various parameter sweeps: shape parameter, speckle size in pixels, pixel crosstalk, and shot noise. In the bottom row, the detrended contrast values are plotted with solid curves and their fitted model (the exponential terms of Eq. (20) only) as black dashed curves. While a parameter is changed, the others remained static as $\alpha = 2$, speckle size $4$ pixels for subplot a) and speckle size $1$ for b-d), shot noise $\sigma _n = 0$, static scatterer ratio $\rho =0$, lateral pixel crosstalk $\sigma _p = 0.5$.

Download Full Size | PDF

4. Correction of overexposure

The goal of this section is to compare the compensation of the contrast distortion of partially saturated images. We demonstrated in the previous sections, that the contrast value decreases as more and more pixels saturate. It can be seen in the detrended contrast curves (Fig. 3), which are the exponential components of the fitted values that the low and high saturation sides can differ significantly ($a_1 \neq a_2$ and $b_1 \neq b_2$). The reason for the difference is that the high and low part of the PDF of the intensity is affected differently by various imaging distortions. This observation suggests that a precise extrapolation cannot be done when near half of the pixels are saturated.

Next, we investigate whether the unsaturated case can be estimated by extrapolation from an artificially generated saturation ratio and contrast curve. Such curves, similar to Fig. 2 and Fig. 3, may be generated from the recorded images by lowering the threshold value step by step, clipping the images and recalculating their contrast maps. The dynamic range of the initial image is limited by the threshold of the sensor as $I_0\in [0, I_{sat}]$. Let us divide this range into $M$ evenly distributed artificial threshold values. After consecutively clipping the original image at each threshold, we get $M$ contrast $\kappa _i$ and saturation ratio $R_i$ maps using Eq. (1), Eq. (18), where $i \in [1, M]$.

We have compared two straightforward function extrapolation strategies. As a reference solution, we fitted the Eq. (20) with $M=7$ threshold levels with constrained nonlinear least-squares optimization. In addition, we explored the near linear behavior of the contrast decreasing in lower contrast cases with a simple and effective linear function. The presented simulation contained the following parameter settings: $\alpha = 1, 2, 5, 8, 20, 50$, speckle size $2$ pixels, shot noise $\sigma _n = 0.01$, lateral pixel cross talk $\sigma _p = 0.5$, static scatterer ratio $\rho =0.2$.

After parameter fitting, each of them was evaluated at $R_s=0$ point to get the unsaturated estimate. The simulation was repeated at 100 illumination levels in the range of $0.25I_{unsat}$ to $5I_{unsat}$. The outcome is presented in Fig. 4(a). The Eq. (20) based extrapolation method resulted in precise, but slow and noisy estimates.

 figure: Fig. 4.

Fig. 4. Comparison of different extrapolation methods: a) contrast and saturation ratio point pairs are fitted by Eq. (20) and evaluated $R_s=0$; b) linear fitting using Eq. (22); c) residual multiplicative error scatter plot and trend of selected, low contrast ($\kappa _s<0.3$) points; d) linear fitting with residual error correction optimized for low contrast values by Eq. (24).

Download Full Size | PDF

Next, the contrast $\kappa _s$ is approximated by a linear function of the saturation ratio and the measured (partially saturated) contrast value:

$$\kappa_{corr, linear} = p_1(\kappa_s, R_s) \cdot R_s + p_0(\kappa_s, R_s).$$

In order to further minimize the computational requirements, we took the advantage that setting the threshold level to 0, the resulting contrast becomes 0 at saturation ratio 1, thus there is no need for further contrast calculations in addition to the already known $\kappa _s$ and $R_s$. Thus this estimation simplifies to:

$$\kappa_{corr, linear} = \frac{\kappa_s}{1-R_s}.$$

The results are shown in Fig. 4(b). As expected, this method produced low estimation error at low contrast values where the contrast function is more linear.

As a continuation, we also analyzed the residual error of the fitting process. Generally, due to the independent distortion sources, error compensation is not feasible, though, for a certain condition range, it is beneficial to reduce the extrapolation imprecision. In order to explore that, multiple simulation setups were used with a wide parameter span. The error was sought as the ratio of the expected and the estimated contrast values $\kappa _{unsat}/\kappa _{corr,linear}$. The relation that describes the error can be derived as the quotient of Eq. (20) and Eq. (22), and by substituting $R_s=0$. The quotient can be approximated with good accuracy by a rational function in the range of $R_s \in [0,1]$ in the following form:

$$\frac{\kappa_{unsat}}{\kappa_{corr,linear}} = \frac{(p_0 + a_1) (1-R_s)}{p_1R_s + p_0 + a_1e^{{-}b_1R_s} - a_2e^{{-}b_2(1-R_s)}} \approx \frac{1+c_1R_s}{1+q_1R_s+q_2R_s^2}.$$

The gathered data and the results of the approximate model (Eq. (23)) is presented in Fig. 4(c) as black dots. For different parameter settings, where the exponential terms of Eq. (20) are small the extrapolation error behaved similarly, independently from the actual contrast value. This observation helps to correct some errors of the simple extrapolation for a wide range of setups. As a rule of thumb, low contrast values ($\kappa _s < 0.3-0.4$, which corresponds to higher flow rate or long exposure time, small speckle size, etc.) satisfy this condition. Including this term in the correction process, the simplified linear extrapolation becomes:

$$\kappa_{corr} = \frac{\kappa_s}{1-R_s} \frac{1+c_1R_s}{1+q_1R_s+q_2R_s^2},$$
where the parameters are $c_1=-0.8, q_1=-0.85, q_2=0.25$ for the selected low contrast region. The corresponding corrected curves can be seen in Fig. 4(d). The simulations of Fig. 3 is repeated using Eq. (24) and can be seen in Fig. 5. The advantage of this form is that it can be evaluated using image operators without pixel-wise fitting routines. Note, to correct the overexposure, $\kappa _s$ and $R_s$ must be calculated for each image, while $c_1, q_1, q_2$ only once. On the other hand, the parameters of the correction term may differ for different setups and need to be customized.

 figure: Fig. 5.

Fig. 5. Linear extrapolation with residual error corrected repetition of the simulations presented in Fig. 3. The dash lines are the expected, unsaturated contrast value levels.

Download Full Size | PDF

5. Measurement

In this section, the application of Eq. (24) is tested in laboratory experiments. A laser spot was generated by an uncollimated single-mode laser diode driven by a constant current source (RLD82PZJ1, 820 $nm$ central wavelength, 220 $mW$, ROHM Co., Ltd., Japan). The laser diode was placed in a mount with a thermoelectric cooling stage (LDM9T, Thorlabs, Newton, NJ, USA) and driven by a current generator (LDP-VRM 01-12 CA, Picolas GmbH, Germany). The images were acquired with a 10 $cm$ focal length and f/11 aperture objective and a monochrome camera of 1536x2048 pixels resolution and $3.45{\times }3.45\;\mathrm{\mu}\textrm {m}^2$ pixel size (Basler ACA2040-55um, Basler Vision Technologies, Germany). Given the linearly polarized laser diode and an additional static and a rotated linear polarizer (Thorlabs, Newton, NJ, USA) the illumination level was varied in wide range.

The first set of measurements observed a static white carton paper and a small area of the backside of a hand with different exposure times (5, 10, 20 $ms$) and changing illumination level. The contrast and Eq. (24) were calculated in a 7x7 pixels sliding window, than averaged for 2x2 $mm^2$ area to make a single measurement dot as shown in Fig. 6. The points shown in the figure belongs to the original and corrected $1/\kappa _s^2$ values. Working with in vivo samples we found that the contrast and saturation ratio curves tends to be linear, with less non-linearity than the ideal scenarios. As Fig. 3 parametric analysis showed, non-idealities separately already result in flattening of the contrast curve. Thus the compensation may be more effective in real applications than the numerical simulations suggest.

 figure: Fig. 6.

Fig. 6. Speckle contrast of the same areas at different illumination levels varied in 1:15 ratio. Red dots scaled on the right axis are the ratio of saturating pixels. Scaled on the left axis, the uncorrected values are shown as blue dots, the corrected values as green dots. Subplots show measurements of 2 $mm$ by 2 $mm$ areas of a) static white paper and b-d) back of a hand with 5, 10, and 20 $ms$ exposure times.

Download Full Size | PDF

As a final experiment, we observed human fingers. 25 frames were taken at $5\:ms$ exposure time first at non-saturating light level, then at three times larger illumination, causing significant area getting overexposured. The spatial contrast is calculated in 11x11 pixels windows and averaged for the frames. Figure 7 shows both the uncorrected and corrected perfusion maps.

 figure: Fig. 7.

Fig. 7. Contrast correction of an overexposured image using Eq. (24). a) the perfusion map generated at medium illumination. The maximum intensity was set to 50% mean value of the sensor’s dynamic range resulting in no saturating pixels. b) perfusion map with overexposured regions using three times higher intensity. Strong distortion can be seen at saturated regions, and better signal to noise ratio at darker areas. c) local mean with color coded saturation (yellow to red areas). d) compensated perfusion map using the proposed linear extrapolation method. The perfusion is estimated as $1/\kappa ^2$.

Download Full Size | PDF

6. Discussion

The computational requirement of the presented linear and corrected extrapolation method is slightly above the standard contrast calculations. This may restrict its real time application in time critical applications, such as to visualize blood flow changes in vivo intraoperative situations [6].

The contrast range reported in practical applications lies in the range of $0.04 < \kappa < 0.6$ restricted by noise, pixel cross-talk, finite laser spatial coherence length and exposition length [31]. Considering these two limiting cases, we showed that up to 50-60% saturating pixel ratio and low contrast range ($\kappa < 0.3-0.4$) simple, linear extrapolation method enables external illumination increase with low blood flow index distortion.

7. Conclusions

Our work suggests that overexposured laser speckle images can be compensated to a wide extent, and common targets can be viewed with better signal-to-noise ratio and more uniform perfusion visualization. Such situations of high dynamic range targets occur during laparoscopic surgeries, various clinical practical situations of diagnosing burn wounds, skin microvasculature, liver, esophagus, and the large intestine. Furthermore, the strict illumination control to avoid any saturation can be waived.

Funding

Eötvös Lóránd Reseach Network (ELKH KO-40/2020); Innovációs és Technológiai Minisztérium (AI National Laboratory Program).

Acknowledgement

Manuscript revision and style correction were done by Ákos Zarándy.

Disclosures

The authors declare no conflicts of interest.

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request. Source code for the presented method is available in Code 1 (Ref. [32]).

References

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

2. P. Li, S. Ni, L. Zhang, S. Zeng, and Q. Luo, “Imaging cerebral blood flow through the intact rat skull with temporal laser speckle imaging,” Opt. Lett. 31(12), 1824–1826 (2006). [CrossRef]  

3. B. Choi, N. M. Kang, and J. Nelson, “Laser speckle imaging for monitoring blood flow dynamics in the in vivo rodent dorsal skin fold model,” Microvasc. Res. 68(2), 143–146 (2004). [CrossRef]  

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

5. W. Heeman, W. Steenbergen, G. M. van Dam, and E. C. Boerma, “Clinical applications of laser speckle contrast imaging: a review,” J. Biomed. Opt. 24(08), 1 (2019). [CrossRef]  

6. A. Mangraviti, F. Volpin, J. Cha, S. I. Cunningham, K. Raje, M. J. Brooke, H. Brem, A. Olivi, J. Huang, B. M. Tyler, and A. Rege, “Intraoperative laser speckle contrast imaging for real-time visualization of cerebral blood flow in cerebrovascular surgery: results from pre-clinical studies,” Sci. Rep. 10(1), 7614 (2020). [CrossRef]  

7. P. H. Carvalho, J. B. Barreto, R. A. Braga, and G. F. Rabelo, “Motility parameters assessment of bovine frozen semen by biospeckle laser (bsl) system,” Biosyst. Eng. 102(1), 31–35 (2009). [CrossRef]  

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

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

10. J. C. Ramirez-San-Juan, C. Regan, B. Coyotl-Ocelotl, and B. Choi, “Spatial versus temporal laser speckle contrast analyses in the presence of static optical scatterers,” J. Biomed. Opt. 19(10), 106009 (2014). [CrossRef]  

11. T. Dragojevic, D. Bronzi, H. M. Varma, C. P. Valdes, C. Castellvi, F. Villa, A. Tosi, C. Justicia, F. Zappa, and T. Durduran, “High-speed multi-exposure laser speckle contrast imaging with a single-photon counting camera,” Biomed. Opt. Express 6(8), 2865–2876 (2015). [CrossRef]  

12. C. Wang, Z. Cao, X. Jin, W. Lin, Y. Zheng, B. Zeng, and M. Xu, “Robust quantitative single-exposure laser speckle imaging with true flow speckle contrast in the temporal and spatial domains,” Biomed. Opt. Express 10(8), 4097–4114 (2019). [CrossRef]  

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

14. C. Liu, K. Kılıç, S. E. Erdener, D. A. Boas, and D. D. Postnov, “Choosing a model for laser speckle contrast imaging,” Biomed. Opt. Express 12(6), 3571–3583 (2021). [CrossRef]  

15. M. Siket, I. Jánoki, K. Demeter, M. Szabó, and P. Földesy, “Time varied illumination laser speckle contrast imaging,” Opt. Lett. 46(4), 713–716 (2021). [CrossRef]  

16. L. Song and D. S. Elson, “Effect of signal intensity and camera quantization on laser speckle contrast analysis,” Biomed. Opt. Express 4(1), 89–104 (2013). [CrossRef]  

17. S. Sunil, S. Zilpelwar, D. A. Boas, and D. D. Postnov, “Guidelines for obtaining an absolute blood flow index with laser speckle contrast imaging,” bioRxiv (2021).

18. J. Zötterman, R. Mirdell, S. Horsten, S. Farnebo, and E. Tesselaar, “Methodological concerns with laser speckle contrast imaging in clinical evaluation of microcirculation,” PLoS One 12(3), e0174703 (2017). [CrossRef]  

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

20. J. Qiu, Y. Li, Q. Huang, Y. Wang, and P. Li, “Correcting speckle contrast at small speckle size to enhance signal to noise ratio for laser speckle contrast imaging,” Opt. Express 21(23), 28902–28913 (2013). [CrossRef]  

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

22. Y. Su, Q. Zhang, and Z. Gao, “Statistical model for speckle pattern optimization,” Opt. Express 25(24), 30259–30275 (2017). [CrossRef]  

23. J. Qiu, P. Li, W. Luo, J. Wang, H. Zhang, and Q. Luo, “Spatiotemporal laser speckle contrast analysis for blood flow imaging with maximized speckle contrast,” J. Biomed. Opt. 15(1), 016003 (2010). [CrossRef]  

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

25. J. W. Goodman, “Statistical properties of laser speckle patterns,” in Laser speckle and related phenomena, (Springer, 1975), pp. 9–75.

26. P. Földesy, M. Siket, I. Jánoki, K. Demeter, and Á. Nagy, “Ensemble averaging laser speckle contrast imaging: statistical model of improvement as function of static scatterers,” Opt. Express 29(18), 29366–29377 (2021). [CrossRef]  

27. S. Covo and A. Elalouf, “A novel single-gamma approximation to the sum of independent gamma variables, and a generalization to infinitely divisible distributions,” Electron. J. Statist. 8(1), 894–926 (2014). [CrossRef]  

28. P. G. Moschopoulos, “The distribution of the sum of independent gamma random variables,” Ann. Inst. Stat. Math. 37(3), 541–544 (1985). [CrossRef]  

29. L. Song, Z. Zhou, X. Wang, X. Zhao, and D. S. Elson, “Simulation of speckle patterns with pre-defined correlation distributions,” Biomed. Opt. Express 7(3), 798–809 (2016). [CrossRef]  

30. H. Spahr, C. Pfäffle, S. Burhan, L. Kutzner, F. Hilge, G. Hüttmann, and D. Hillmann, “Phase-sensitive interferometry of decorrelated speckle patterns,” Sci. Rep. 9(1), 11748 (2019). [CrossRef]  

31. D. D. Postnov, X. Cheng, S. E. Erdener, and D. A. Boas, “Choosing a laser for laser speckle contrast imaging,” Sci. Rep. 9(1), 2542 (2019). [CrossRef]  

32. P. Földesy, M. Siket, I. Jánoki, and Á. Nagy, “Correction of overexposure in laser speckle contrast imaging,” GitHub (2022), https://github.com/fpeter71/OverExposureLSCICompensation.

Supplementary Material (1)

NameDescription
Code 1       The source code for the presented method

Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request. Source code for the presented method is available in Code 1 (Ref. [32]).

32. P. Földesy, M. Siket, I. Jánoki, and Á. Nagy, “Correction of overexposure in laser speckle contrast imaging,” GitHub (2022), https://github.com/fpeter71/OverExposureLSCICompensation.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (7)

Fig. 1.
Fig. 1. Density functions of 256x256 pixels of simulated speckle patterns are shown. The used parameters are $\alpha = 4$ , speckle size 4 pixels, shot noise $\sigma _n = 0.01$ , lateral pixel crosstalk $\sigma _p = 0.5$ , static scatterer $\rho =0$ . Subplot b) shows the speckle pattern at maximum light intensity ( $I_{unsat}$ ) where no saturation occurred. The subplots a), c), d) correspond to $0.25, 2.5, 4$ times this intensity. The contrast and saturation ratio of the four speckle patterns are $K_s = 0.46, 0.47, 0.41, 0.27$ and $R_s = 0\%, 0\%, 9.1\%, 43.5\%$ , respectively.
Fig. 2.
Fig. 2. Subplot a) shows the contrast values and the saturation ratio as a function of the mean observed illumination. The yellow markers corresponds to the numerical simulation of Fig. 1. Subplots b) and c) show the linear and the exponential components of the contrast approximation as a function of the saturation ratio.
Fig. 3.
Fig. 3. The upper subplots show the contrast value as a function of the saturation ratio under various parameter sweeps: shape parameter, speckle size in pixels, pixel crosstalk, and shot noise. In the bottom row, the detrended contrast values are plotted with solid curves and their fitted model (the exponential terms of Eq. (20) only) as black dashed curves. While a parameter is changed, the others remained static as $\alpha = 2$ , speckle size $4$ pixels for subplot a) and speckle size $1$ for b-d), shot noise $\sigma _n = 0$ , static scatterer ratio $\rho =0$ , lateral pixel crosstalk $\sigma _p = 0.5$ .
Fig. 4.
Fig. 4. Comparison of different extrapolation methods: a) contrast and saturation ratio point pairs are fitted by Eq. (20) and evaluated $R_s=0$ ; b) linear fitting using Eq. (22); c) residual multiplicative error scatter plot and trend of selected, low contrast ( $\kappa _s<0.3$ ) points; d) linear fitting with residual error correction optimized for low contrast values by Eq. (24).
Fig. 5.
Fig. 5. Linear extrapolation with residual error corrected repetition of the simulations presented in Fig. 3. The dash lines are the expected, unsaturated contrast value levels.
Fig. 6.
Fig. 6. Speckle contrast of the same areas at different illumination levels varied in 1:15 ratio. Red dots scaled on the right axis are the ratio of saturating pixels. Scaled on the left axis, the uncorrected values are shown as blue dots, the corrected values as green dots. Subplots show measurements of 2 $mm$ by 2 $mm$ areas of a) static white paper and b-d) back of a hand with 5, 10, and 20 $ms$ exposure times.
Fig. 7.
Fig. 7. Contrast correction of an overexposured image using Eq. (24). a) the perfusion map generated at medium illumination. The maximum intensity was set to 50% mean value of the sensor’s dynamic range resulting in no saturating pixels. b) perfusion map with overexposured regions using three times higher intensity. Strong distortion can be seen at saturated regions, and better signal to noise ratio at darker areas. c) local mean with color coded saturation (yellow to red areas). d) compensated perfusion map using the proposed linear extrapolation method. The perfusion is estimated as $1/\kappa ^2$ .

Equations (24)

Equations on this page are rendered with MathJax. Learn more.

κ s 2 = σ 2 ( I s ) I s 2 ,
υ 2 ( T , τ c ) = β T 2 0 T 0 T ρ | g 1 n ( t t ) | 2 + ( 1 ρ ) | g 1 n ( t t ) | 2 d t d t ,
I s = ( 1 ρ ) S ( λ S ) + ρ D ( α D , β D ) + W ,
f ( I i ) = β α Γ ( α ) I i α 1 e β I i ,
F ( I i ) = γ ( α , β I i ) Γ ( α ) ,
F s a t ( I i ) = { F ( I i ) , if I i < I s a t 1 , otherwise . ,
E ( γ ( α , β I s a t ) ) = 1 β γ ( α + 1 , β I s a t ) γ ( α , β I s a t )
E ( f s a t ) = α β γ ( α + 1 , β I s a t ) Γ ( α + 1 ) + I s a t ( 1 F s a t ( I i ) )
σ 2 ( f s a t ) = α β 2 γ ( α + 2 , β I s a t ) Γ ( α + 1 ) E ( f s a t ) 2 + I s a t 2 ( 1 F s a t ( I i ) ) .
κ 2 κ s 2 ( 1 + I s a t E ( f ) ( 2 b 1 ) ) .
R = 1 F s a t ( I i )
Y = E ( f s a t ) I s a t R = α β γ ( α + 1 , β I s a t ) Γ ( α + 1 )
Let       A = 2 Y I s a t 2 = 1 Γ ( α ) β α e β I s a t [ α I s a t α 1 I s a t α β ] ,
and let       B = R I s a t = 1 Γ ( α ) β α I s a t α 1 e β I s a t ,
then       A B = α + β I s a t ,
I s a t A B = β ,
( I s a t A B ) I s a t A B = α
R s = b i n ( I i ) , where
b i n ( I i ) = { 0 , if I i < I s a t 1 , otherwise .
κ s ( R s ) p 1 R s + p 0 + a 1 e b 1 R s a 2 e b 2 ( 1 R s ) .
κ c o r r , l i n e a r = p 1 ( κ s , R s ) R s + p 0 ( κ s , R s ) .
κ c o r r , l i n e a r = κ s 1 R s .
κ u n s a t κ c o r r , l i n e a r = ( p 0 + a 1 ) ( 1 R s ) p 1 R s + p 0 + a 1 e b 1 R s a 2 e b 2 ( 1 R s ) 1 + c 1 R s 1 + q 1 R s + q 2 R s 2 .
κ c o r r = κ s 1 R s 1 + c 1 R s 1 + q 1 R s + q 2 R s 2 ,
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.