Robust quantitative single-exposure laser speckle imaging with true flow speckle contrast in the temporal and spatial domains

A systematic and robust laser speckle contrast imaging (LSCI) method and procedure is presented, covering the LSCI system calibration, static scattering removal, and measurement noise estimation and correction to obtain a true flow speckle contrast and the flow speed from single-exposure LSCI measurements. We advocate to use as the speckle contrast instead of the conventional contrast K as the former relates simply to the flow velocity and is with additive noise alone. We demonstrate the efficacy of the proposed true flow speckle contrast by imaging phantom flow at varying speeds, showing that (1) the proposed recipe greatly enhances the linear sensitivity of the flow index (inverse decorrelation time) and the linearity covers the full span of flow speeds from 0 mm/s to 40 mm/s; and (2) the true flow speed can be recovered regardless of the overlying static scattering layers and the type of speckle statistics (temporal or spatial). The fundamental difference between the apparent temporal and spatial speckle contrasts is further revealed. The flow index recovered in the spatial domain is much more susceptible to static scattering and exhibit a shorter linearity range than that obtained in the temporal domain. The proposed LSCI analysis framework paves the way to estimate the true flow speed in the wide array of laser speckle contrast imaging applications.


Introduction
When coherent light interacts with a turbid medium, the interference between the outgoing waves produces grainy speckle patterns which encode the phase fluctuation of all rays (random phasors) reaching a point.The contrast of laser speckles reduces with the motion of scatterers inside the turbid medium.Laser speckle contrast hence can be used to infer the dynamic property of the medium.Laser speckle contrast imaging (LSCI, see recent reviews [1][2][3][4]) has now been widely used in monitoring blood flow in brain, skin, retina, arthrosis and etc due to advantages including simplicity, high spatial and temporal resolution, and large field of view without scanning [5][6][7][8].
Although LSCI has a wide range of applications and a long history, the recovery of absolute flow velocity from LSCI measurements remains a challenge, especially when the measurement is compounded by static scattering and noise.For static scattering in laser speckle imaging, Li et al. showed that the static scattering effect can be partially suppressed by using the temporal rather than spatial contrast analysis of laser speckles [9] as the static scattering is an invariant quantity with time.Zakhraov et al. [10,11] presented a data processing scheme to correctly separate dynamic and static components within the speckle contrast based on their different decorrelation behaviour across speckle patterns captured at consecutive times.Dunn et al. [6,[12][13][14] later demonstrated a multi-exposure laser speckle contrast imaging method, which quantifies and eliminates the influence of static scattering from speckle contrasts measured under different exposure times using a laser speckle contrast model.For LSCI measurement noise, the correction of the variance of the shot noise and sensor dark currents were found to be crucial to estimate the true speckle contrast [15,16].Yuan et al. [16] increased the signalto-noise ratio (SNR) of LSCI with noise correction to detect small blood flow changes caused by brain activity.A systematic study and recommended practical recipe to obtain true flow velocity from LSCI measurements addressing both static scattering and measurement noise is, however, still lacking.
In this article, we analysed laser speckle flow imaging from the first principle and provided a complete procedure covering the LSCI system calibration, static scattering removal, and measurement noise estimation and correction to obtain a genuine flow speckle contrast and the flow speed from single-exposure LSCI measurements.We demonstrated the power of our LSCI analysis recipe by imaging phantom flow at varying speeds.Experimental results show that our procedure greatly enhances the linear sensitivity of the flow index (defined as the inverse decorrelation time) and the linearity covers the full span of flow speeds from 0 mm/s to 40 mm/s.The true flow speed is recovered regardless of the overlying static scattering layers and the type of statistics (temporal or spatial).The proposed LSCI analysis framework hence paves the way to estimate the true flow speed in the wide array of laser speckle contrast imaging applications.

Theoretical basis
The spatial intensity distribution of the speckle pattern fluctuates with the motion of the scattering particles under the illumination of coherent light.The recorded pattern by a camera is the integration of all instantaneous speckles over the exposure time.The faster the scattering particles move, the more blurred the recorded pattern becomes.The degree of blurring is quantified by the contrast [17] given by 2 () where I represents the mean of light intensity I over a small region (spatial contrast) or over a short durance of time (temporal contrast).For "fully developed" static speckles, the spatial contrast K equals to 1.
We will assume the scattered electric field containing both dynamic and static components tt

( ) [ ( ) ]
ii fs with ω being the angular frequency of light.The dynamic component consists of photons which have at least been scattered by moving scatterers (flow) once and the static component consists of photons being scattered by static scatterers alone.The electric field temporal autocorrelation function can be written as where means average over t,

 
is the electric field temporal autocorrelation function related to flow, and ( ) ( ) / (0) g GG

 
, the full Siegert relation [18,19] expresses where ()   is a parameter that accounts for the reduction in the measured contrast due to averaging (by the detector) over uncorrelated speckles.Note 1 g is real and non-negative where the average intensity   .We will use 2 K as the speckle contrast instead of the conventional contrast K as the former relates simply to the flow velocity and is with additive noise alone.
Equation ( 5) is the theoretical temporal contrast from the random process taken by the photons migrating through a turbid medium.Some complexities arise when evaluating 2  K from measurement data.First the measurement noise (of zero mean) introduces an extra variance term 2 noise  .Second when using spatial ensemble average for the evaluation of the contrast rather than temporal statistics, the extra terms appearing in 2  K will be  43 where / c xT   and c  is the decorrelation time.The inverse decorrelation time increases with the flow speed and can be regarded as the flow index.
In the next section, we will examine system calibration, sample measurement and data analysis to provide a complete procedure for static scattering removal, and measurement noise estimation and correction to obtain a true flow speckle contrast and the flow speed from singleexposure LSCI measurements.

System calibration, measurement, and data analysis
Let's consider a set of speckle images i I ( Here the recorded image consists of the static component s I , the dynamic component f I and the noise n.The noise [15] presented in the measurement is mainly comprised of the dark counts The dark counts d n and the variance of dark counts 2 2 () Var n n n  can be easily acquired by taking multiple dark frames at the same exposure time and camera gain (with all light off).One could use the temporal average to get the dark count and its variance pixel by pixel when the number of the dark frames is large ( 50  ) or use the spatial average over a sliding PP NN  (typical P N =7) pixel window otherwise.If the behavior of dark counts is assumed uniform across the whole sensor frame, the mean and the variance of the dark counts are given by further averaging over the whole sensor frame.In many cameras, the recorded intensity has been pre-subtracted by certain base.In this case, d n replaced by zeros.We would always subtract the dark counts from i I before further analysis.After this pre- processing, the speckle image becomes ( , ) ( , ) ( , ) . The noise satisfies 0 Var n Var n Var n .
The shot noise obeys a Poisson distribution with the mean 0 si n  and a variance equal to ( , ) / c i I x y  as the camera converts the photoelectrons to digital counts where  is the analog to the digital conversion factor [20]  .The  factor is typically the same across the sensor frame and thus is obtained by further averaging over the sensor frame.Under such a shot noise model, where an extra quadratic term in ( , )  where the coherence factor of the imaging system is defined as

System calibration
As stated earlier, the dark counts d n and the variance of dark counts () Var n is first acquired by taking multiple dark frames at the identical experimental condition (the same exposure time, camera gain and etc) with all light off.The other system parameters (the coherence factor  and the behavior of   c Var n vs c I ) can be directly evaluated from a set of fully developed speckle images taken on a pure static scattering sample such as a reflection standard.Indeed, according to Eq. ( 17 for 0  in this case.Here the spatial average should use the largest window size (full image if possible) due to the reason discussed in [21] for the temporal speckle contrast.For the spatial speckle contrast P N in Eq. ( 19) should take the same value used for the contrast calculation.
The  factor is typically the same across the sensor frame and thus obtained by further averaging over the sensor frame.
The noise variance

  c
Var n associated with a particular c I can be found as

Sample measurement and data processing
The sample containing both dynamic and static scatterers are then imaged under the identical experimental condition to yield a new set of dark current removed speckle images ( , ) The dynamic fraction  can be determined using Eq. (17)as for The noise level can also be estimated directly from the set of speckle images via in the spatial domain (they are equal by ergodicity).Using Eq. ( 13) and ( 16), the temporal speckle contrast is then and the spatial speckle contrast is then from which we can identify the non-ergodic contribution from the static field to be Finally, the velocity information of the sample can be obtained by solving the flow contrast 2 f K and fitting to Eq. ( 8) to obtain the decorrelation time c  and the flow index.
Figure 1 outlines the complete procedure of system calibration, sample measurement, noise correction, and static scattering removal for robust quantitative single-exposure laser speckle imaging.System calibration first determines the camera dark current , the coherent factor  , and the noise parameters  and  from measuring fully developed speckles produced by a pure static reflection standard.The true flow contrast 2 f K of dynamic samples are afterwards obtained by removing the static scattering and correcting the measurement noise from the measured temporal or spatial speckle contrasts.The flow decorrelation time and speed are then be determined from 2 f K with a proper flow velocity model such as Eq. ( 8).   Figure 3 shows the dark current of the camera with a distribution centered at 0. This means that the dark current of the camera has been pre-subtracted and d n should be set to 0.

Results of system calibration
A set of 150 reflectance images from the reflectance standard were then recorded.The coherence factor  of the system was then computed with Eq. ( 19) for different step size  (see Fig. 4).The coherence factor  reduces slightly with  owing to the inevitable system instability.The proper system coherence factor is obtained by extrapolating to =0.3  (=0.540/67, determined by the exposure time 40 msec and the acquisition time 67 msec), yielding  =0.3144.
Fig. 4 The coherence factor  reduces with the step size  .
By varying the intensity attenuator and the reflection ratio of DMD, multiple sets of 150 reflectance images from the reflectance standard were recorded.The noise variance was computed with Eq. ( 20) or (21) in the temporal or spatial domain.The noise variance computed with either approach agrees with each other.The noise variance in the spatial domain, however, has lower standard error and is preferred (see Fig. 5).The computed noise variance increases with the step size  and the light intensity.The proper noise variance is obtained by extrapolating to =0.3  as well as the determination of  above.Figure 6 shows the noise characteristics of the imaging system by extrapolating to =0.3  .The shadow represents the error range given by the standard deviation computed from five separate sets of measurements.Table 1 displays the noise parameters of () c i Var n by fitting with Eq. (10).was often used.The value of this factor calculated from the measurement is observed to decrease with the light intensity (see Fig. 7).The assumption of a constant '  is thus not correct, attributed to the nonzero  mainly caused by the light source fluctuations.Fig. 7 The analog-to-digital factor '  decreases with the light intensity.

Importance of static scattering removal and noise correction
A stack of 250 images were taken for the flow phantom at each flow speed ranging between 0 and 40 mm/s.The dynamic fraction  was computed with Eq. ( 22).Fig. 8(a) shows the 2D distribution of  with an average value of 0.871 over a region of interest (ROI) when flow speed is 0 (Brownian motion alone).The extracted value of  stays unchanged when the flow speed increases (see Fig. 8(b)).The non-uniformity of the dynamic fraction is caused by the imperfect glass tube.
Fig. 8 The dynamic fraction  over a ROI. (a) 2D distribution when the flow speed is 0. (b) The average dynamic fraction versus the flow speed.
The noise contrast 2 noise  was computed using Eq. ( 23) or Eq. ( 24) in the temporal or spatial domain, respectively.Fig. 9 shows the temporal and spatial 2 noise  when the flow speed is 0. The temporal 2 noise  has much higher spatial resolution than the spatial one.The average temporal    The flow speed is directly related to the decorrelation time.The inverse decorrelation time 1/ c  increases with the flow speed and may serve as its proxy.The inverse decorrelation time extracted from fitting Eq. ( 8) is shown in Fig. 11.The sensitivity of uncorrected 1/ c  to the flow speed is very poor and loses linearity around 5 mm/s whereas the corrected 1/ c  shows excellent linearity over the whole range up to 40 mm/s.The corrected one with both noise and static scattering removal also exhibits the least standard deviation.

Effects of different static scattering
The efficacy of the static scattering removal is then investigated.One part of the glass tube was coated with a scattering layer (dried colloidal suspension) and the same set of the measurements were performed.The region A (=0.83, average light intensity = 670) in the green rectangle is covered by the static scattering layer whereas the region B (=0.88, average light intensity = 810) in the red rectangle is directly exposed (see Fig. 12).Both regions should have identical flow speed.13 (a,d)).After noise correction, the agreement between A and B significantly improves although the discrepancy between their recovered 1/ c  is appreciable (see Fig. 13 (b,e)).With a further static scattering removal, the gap between 1/ c  for the two regions in (e) almost disappeared (see Fig. 13 (c,f)).The degrade in the performance for faster flow speeds is caused by the loss of SNR at higher speeds.The above results show that different static scattering can be successfully removed to obtain the true flow velocities.
Fig. 13 (a-c) The temporal speckle contrast (uncorrected, after noise correction, after both noise correction and static speckle removal) for ROI A and ROI B; (d-f) the recovered corresponding inverse decorrelation time.
To further show the agreement of the flow speed in ROI A and B, the error in the recovered 1/ c  can be directly estimated using the uncertainty in the noise variance.The noise contrast depends on light intensity alone when the imaging system has been specified.At higher speeds, the uncertainty in the noise variance starts to dominate as the flow contrast steadily reduces.Fig. 14 shows the flow speeds in regions A and B indeed agree with each other within the system uncertainty given in Fig. 6 and Table 1.
Fig. 14 The flow speed at region A and B agrees with each other within system uncerntainty.

Temporal speckles vs spatial speckles
The speckle contrast analysis can not only be performed within the temporal domain presented in Sec.3.3.1 and 3.3.2but also in the spatial domain.The two different approaches have their own merits.
Fig. 15 (a-c) The temporal and spatial speckle contrast (uncorrected, after noise correction, after both noise correction and static speckle removal) for ROI A; (d-f) the recovered corresponding inverse decorrelation time.
Figure 15 compares the temporal and spatial speckle contrast and the inverse decorrelation time for ROI A. The uncorrected temporal and spatial speckle contrasts and inverse decorrelation times differ significantly caused by the non-ergodic static scattering 2 ne  (see Fig. 14 (a,d)).After noise correction, the temporal speckle contrast and inverse decorrelation time performs much better than the spatial counterparts which still retain 2 ne  (see Fig. 14 (b,e)).
With a further static scattering removal, 2 f K and 1/ c  in the temporal and spatial domains tend to agree with each other (see Fig. 13 (c,f)).However, a careful examination of the recovered 1/ c  in (f) reveals their difference.The inverse decorrelation time recovered in the spatial domain shows much less variation yet the linearity range of the temporally recovered 1/ c  expands to much higher speeds.The latter behavior can be attributed to the difficulty of static scattering removal inside the spatial domain where a subtraction between the measured contrast

2D flow profile
In addition to the above LSCI analysis of the overall behavior of the flow contrast and inverse decorrelation time versus flow speed, the result of LSCI imaging of a specific 2D region is shown in Fig. 16 (v=3 mm/s).The flow speed is observed to increase closer to the center of the tube.The flow speed cross sectional profile 1/ c  marked in Fig. 16 (a) fits well by a Newtonian flow profile [22].Comparing speckle contrasts in the temporal and the spatial domain, the latter contains one additional term of 2 ne  .This fact explains the apparent increase of the spatial speckle contrast with the static scattering [23].The inverse decorrelation time recovered in the spatial domain shows much less variation yet with a much shorter linearity range than that obtained in the temporal domain.The rapid deterioration of the performance of 2 f K in the spatial domain is caused by the difficulty of static scattering removal which requires a subtraction between the measured spatial speckle contrast 2  K and 2 ne  .At higher speeds, the error in 2 ne  dominates and the flow speckle contrast 2 f K can no longer be accurately estimated.This observation is fundamental in selecting the appropriate statistics in analysing the LSCI measurements.A general guideline is that the spatial speckle contrast should be avoided when 2 f K is not larger than the error in 2 ne  .
A Lorentz velocity distribution model for the moving particles is assumed to relate 2 f K to the flow speed in our study.Different velocity distribution models may be assumed [24].However, in the typical situation of much longer exposure time compared to the decorrelation time (as in our study), the relation Eq. ( 8) stills holds other than a trivial pre-factor.The  is the vacuum wavelength of the incident light and n is the refractive index of the medium.For example, it yields 2.0 0.1  mm/s when the input flow speed is 5 mm/s and 4.2 0.2  mm/s when the input flow speed is 10 mm/s (see Fig. 15).
Finally, although our study is on the single-exposure laser speckle imaging, the same analysis methodology can be carried over to the multiple-exposure LSCI.The latter gains the advantage over the former when probing the flow velocity distribution.However, the much simpler and faster single-exposure LSCI performs as well as multiple-exposure LSCI when the detailed velocity distribution is not of interest as long as the exposure time is much larger than the decorrelation time and the outlined analysis procedure is followed.
The code for the proposed analysis procedure has been provided in GitHub [25].

Conclusion
In summary, a systematic and robust laser speckle flow imaging method and procedure has been presented, covering the LSCI system calibration, static scattering removal, and measurement noise estimation and correction to obtain a true flow speckle contrast and the flow speed from single-exposure LSCI measurements.The power of our LSCI analysis recipe has been demonstrated by imaging phantom flow at varying speeds, showing that (1) our recipe greatly enhances the linear sensitivity of the flow index and the linearity covers the full span of flow speeds from 0 to 40 mm/s; and (2) the true flow speed is recovered regardless of the overlying static scattering layers and the type of speckle statistics (temporal or spatial).The difference and merits of the temporal and spatial speckle contrasts have been compared and a guideline for selecting the appropriate statistics for LSCI has been provided.The proposed LSCI analysis framework paves the way to estimate the true flow speed in the wide array of laser speckle contrast imaging applications.

dn
and the shot noise s averaging over the sliding window of size P N .Equation (17) is also correct for a pure static scattering sample producing fully developed speckles (ρ = 0) as long as ji  .
or spatial statistics, respectively.Multiple sets of such speckle images under the identical experimental condition and varying incident intensities are measured to cover the full range of c I .By fitting to Eq. (10), the system noise behavior is then characterized.We note the above results on  and dependence on ∆ can be observed due to the inevitable system instability.In this case, the correct values of  and

Fig. 1
Fig.1Experimental and data analysis framework for robust quantitative single-exposure laser speckle imaging.

Figure 2
Figure 2 shows the schematic diagram of the experimental setup.Light from a DPSS red laser (LSR671ML, λ = 671nm, Lasever, Ningbo, China) illuminated the sample and the speckle images were recorded by a 12bit camera (MER-125-30UM, Daheng Imaging, China, 1292×964 pixels, 3.75μm×3.75μm)with an exposure time set between 20 and 40 msec.The DMD (DLC9500P24 0.95VIS) acted as a reflection mirror here.In system calibration, light reflectance from a Lambertian reflection standard was recorded with the exposure time set at 40 msec and a total of 150 images captured at a frame rate of 15 fps.The system characteristics under different levels of light illumination was obtained by varying the intensity attenuator and the reflection ratio of DMD.In flow velocity measurement experiments, Intralipid-2% suspension (scattering coefficient = 1.7 1 mm  ) inside a glass tube (inner diameter 1mm, outer diameter 2mm) is used to simulate blood flow.The flow rate in the glass tube is set by adjusting the driving speed of the fuel injection pump, covering the whole range from 0 to 40 mm/s in this study.A stack of 250 raw speckle images of the dynamic sample was acquired with an exposure time set at 40 msec and a frame rate of 15 fps for each flow speed.

Fig. 2
Fig. 2 Schematic of the experimental setup.

Fig. 3
Fig. 3 Dark current of the camera.

Fig. 5
Fig. 5 Noise variance increases with the step size  and the light intensity for (a) two particular intensities, and (b) ∆=0.3 and 1.

Fig. 6 3 .
Fig. 6 Noise characteristics of the imaging system: (a) Noise variance and (b) noise contrast 2 noise  versus previous LSCI experiments, an analog-to-digital conversion factor[9]

Figure 10
Figure10shows the temporal speckle contrast2  K computed from the data set (original, after noise correction, after both noise correction and static scattering removal yielding 2 f K ).

Fig. 10
Fig. 10 Correction of the temporal speckle contrast 2 K . 2 K after both noise correction and static scattering removal yields 2 f K .

Fig. 11
Fig. 11 The inverse decorrelation time from (a) uncorrected and (b) corrected temporal speckle contrast.

Fig. 12
Fig. 12 ROI A and B (covered with an extra static scattering layer) are imaged.

Figure 13
Figure 13 compares the temporal speckle contrast and the inverse decorrelation time for ROI A and B. The inverse decorrelation time from the uncorrected speckle contrast differs

2 ne dominates and the flow speckle contrast 2 fK
is required.At higher speeds, the error in computed in the spatial domain fails to obtain the true flow speed.

Fig. 16
Fig. 16 (a) ROI selected for analysis.(b) Flow index 1/ c  for the ROI. (c) The flow speed cross sectional

2 fK
The measured temporal and spatial speckle contrasts for flow imaging are affected by both the presence of static scattering and measurement noise.Their values always differ from each other except for a pure dynamic medium without static scattering.The spatial speckle contrast contains one extra term 2 ne  due to the non-ergodic static scattering than the temporal counterpart.Nevertheless, a common true flow speckle contrast can be defined in both the temporal and the spatial domains.A complete procedure covering the LSCI system calibration, static scattering removal, and measurement noise estimation and correction to obtain the true flow speckle contrast 2 f K and the flow speed from single-exposure LSCI measurements has been detailed here.The recovered inverse decorrelation time 1/ c  from 2 f K exhibits excellent linearity against the flow speed over the full span from 0 to 40 mm/s.The true 1/ c  is obtained regardless of the overlying static scattering layers and the type of measured contrasts (temporal or spatial speckle contrasts).
compares favorably with the input value where

Table 1
Fitted noise parameters of