1. Introduction

Laser Speckle Contrast Imaging (LSCI) is a versatile, low-cost modality for measuring blood flow in exposed tissue, commonly applied to measuring cerebral blood flow (CBF) [1–7]. It is particularly attractive as an intraoperative tool for measuring CBF, as it provides full-field relative blood flow speed maps with the use of only a coherent light source and a camera, both easily adaptable with typical microscopes used in neurosurgery [8]. Since LSCI is based on the extraction of quantitative information from the local spatial speckle contrast in an image, K = σ_I/〈I〉 [9], it is sensitive to image blurring from externally induced movements. Consequently, LSCI is vulnerable to CBF overestimation from brain tissue movements. These movement-related artefacts may be due to the pulsatile nature of blood flow, breathing-related movement, longer-timescale transients due to tissue swelling and deformation, and non-physiological artefacts due to mechanical vibration or motion of the imaging systems components. The aforementioned artefacts can be classified into two categories: those occurring at a speed comparable to the integration time of the camera and those that take place on a longer timescale. The former category manifests itself as spikes in the extracted flow speeds and therefore does not contain physiologically relevant information. The work in reference [8] shows a method for discarding such non-physiological flow information. The latter category can be further classified into two types of artefacts that can be corrected for either in real-time or in post-processing: lateral movement and defocusing. Lateral movement can be corrected using a variety of image registration algorithms in post-processing [10,11]. Correcting for defocusing is more challenging since post-processing algorithms cannot bring the image back into focus reliably [12]. A variety of autofocus algorithms have been reported [13–17]; these algorithms strive to maximize contrast in either the intensity domain or the spatial frequency domain. Applying such methods to either the un-processed speckle images or LSCI contrast maps is challenging. Indeed, applying any traditional algorithms to the un-processed speckle images is not feasible because a) near-IR (NIR) light that is typically used for LSCI measurements, due to its low absorption in tissue [18], shows little intrinsic contrast in the un-processed images due to the very similar optical absorption by the tissue and blood vessels, b) the abundance of high spatial frequency components in the speckle images makes spatial frequency domain analysis not feasible. Applying the autofocus algorithms to the LSCI processed image directly, e.g. by maximizing the measured speed in the vessel, has two key drawbacks. First, the LSCI image is inherently noisy (both spatially and temporally), and typically requires a moving average filter for reliable extraction of quantitative data. Secondly, intrinsic changes in the signal from neurovascular coupling, falling blood pressure, or an ischemic insult, may result in a false trigger to autofocus, thus making the system less robust.

As such, in this paper we propose a different measure of misfocus (i.e., axial deviation from the plane of best focus) - the kurtosis of the vessel flow profile - as a more sensitive and reliable measure of defocusing compared to direct measurements of the LSCI signal magnitude (i.e., relative flow speed values). We show that the kurtosis of the flow profile is a more sensitive, less noisy measure of the imaging-system focal-plane location relative to an imaged vessel (i.e., misfocus) and demonstrate a system and method for using our proposed measure to correct for defocusing in real-time in a controlled experiment, without using a moving average temporal filter. This paper is organized as follows: section 2 (Methodology) describes rationale and background for the measure we used; section 3 (Results) presents our experimental setup and the results of our characterization experiments; and section 4 (Discussion) addresses practical concerns in applying our proposed measure and system to in-vivo intra-operative and chronic optical recordings.

2. Methodology

This section describes the principle of operation and computational algorithm in our proposed measure for misfocus.

2.1. Kurtosis of flow profile

Kurtosis is formally defined as the fourth standardized moment of a particular distribution [19]:

In the above equation, X is a distribution, μ is the mean, and E is the expected value. This measure describes the deviation of the distribution from a normal one, alternatively describing distribution’s “peakedness”. While Eq. (1) is a statistical measure that is a property of a distribution, we found it to be a good descriptor of the behaviour of vessel flow profiles as they become defocused. Qualitatively, defocusing acts as a spatial low-pass filter that acts on the LSCI-derived flow speed (speckle flow index, SFI) maps, causing a blur of high spatial frequency components that are predominantly found on the edges of the vessel. As a result, the flow profile of a defocused vessel looks more “peaky”, as shown in Fig. 1. By mathematically treating the transverse SFI profile, v(x), where x is the coordinate in the transverse direction to the flow, as a distribution, we compute its kurtosis as:

 

Fig. 1 (a) Vessel kurtosis as a function of misfocus (vertical position, z, of the objective compared to the plane of best focus) for an arbitrary vessel in a mouse brain. Images were taken using an Olympus microscope (model BX61WI) using a 4× objective, NA = 0.1. A 680 nm vertical-cavity surface-emitting laser was used for illumination, while the camera exposure time was set to 5 ms. (b) (i) - (vii) show the normalized relative flow speed maps of the vessel for the different misfocus values indicated in (a), (x is the transverse direction across the vessel). (c) (i) - (vii) show the normalized relative flow speed profiles obtained by averaging the results in (b) along the columns.

Download Full Size | PDF

The resulting curve of kurtosis as a function of misfocus has the form shown in Fig. 1. As indicated in Fig. 1, the ability of Eq. (2) to accurately correlate with misfocus magnitude breaks down when the vessel profile is visibly defocused and the profile can no longer be accurately viewed as a distribution. As such, we augment the measure of misfucus by multiplying it by the “apparent width” of the vessel, which we quantify in the next subsection. We proposed that this combined measure, which we term ζ, serves as a more accurate description of misfocus than evaluating the integrated SFI signal directly.

2.2. Computational routine

Figure 2 illustrates the process of computation of ζ. It is inherent from the form of Eq. (2) that computing kurtosis requires a selection of a vessel within the region of interest (ROI) (Fig. 2). To estimate the flow profile, a vessel is first cropped (Fig. 2(i)), and placed within a circular field (Fig. 2(ii)). The vessel’s orientation was calculated by computing the covariance matrix, Σ, of the cropped image (by assuming a bi-variant distribution of pixel values as data points) (Fig. 2(iii)). The eigenvector associated with the smaller eigenvalue corresponded to the vessel’s cross-sectional profile direction (cross-variance), while the eigenvector associated with the larger eigenvalue corresponded to the longitudinal direction along the vessel. The flow map of the vessel was then rotated using bilinear interpolation (Fig. 2(iv)), and the cross-sectional flow profile from which the kurtosis was computed was obtained by averaging along pixel columns (Fig. 2(v)). Lateral movements were corrected by sequential cropping of the flow profile around the profile mean such that the flow profile sampling region was approximately centred on the vessel (Fig. 2(v), dashed red lines). This was necessary due to the inherent noise of fourth order statistics which would produce over-estimations of kurtosis for lateral shifts. Kurtosis of the profile was computed (Fig. 2(vi)), and then multiplied by the negative of the magnitude of the smaller eigenvalue of Σ to produce ζ (Fig. 2(vii)), our chosen measure of misfocus. This smaller eigenvalue is indicative of the pixel intensity distribution in the direction transverse to the vessel, and therefore can be seen as a measure of the “apparent width” of the vessel. The above operations were performed using an efficient MATLAB routine which did not limit the speed of our image acquisition or correction for misfocus.

 

Fig. 2 Algorithm for computing misfocus measure ζ: initially, a square ROI is chosen for a given vessel in the field of view of the SFI map. Then, for each frame (i) the square ROI is cropped, (ii) the ROI is multiplied by a circular mask to improve the calculation of vessel orientation, (iii) a covariance matrix, Σ, is computed for the circular ROI and then Σ is eigen-decomposed, (iv) the eigenvectors of Σ, which reflect vessel orientation, are used to rotate the image in (ii) via a bilinear interpolation, (v) the cross-sectional flow profile is extracted from averaging along the pixel columns. The profile is cropped and re-centred to normalize-out small in -plane displacements (vi), the kurtosis is computed as per Eq. (2) for the given depth (red circle), (vii) the kurtosis valueis multiplied by the negative of the smaller Σ (i.e., the cross-variance of circular ROI – the “apparent width” of the vessel) to arrive at ζ. The maximum value of ζ corresponds to the location of best focus for the vessel.

Download Full Size | PDF

The measure ζ was chosen as it was observed to consistently decrease with the distance of the vessel from best focus for every vessel analysed. In contrast, the profile kurtosis increases with the distance of the vessel from best focus until the vessel image becomes blurry, beyond which the profile kurtosis decreases with increasing misfocus. This limits the axial range in which kurtosis is a robust measure of misfocus. The cross-variance, on the other hand, has an observed plateau-like feature when the vessel is near best focus. Combining these two measures into our measure ζ avoids this limitation without increasing measurement noise.

Since ζ is a product of two functions: kurtosis, which has three observed extrema, and cross-variance, which has one observed extremum, it is possible, in principle, for ζ to have local maxima in addition to a global maximum. In practice, however, no local maxima are observed. Qualitatively, this is because significant change in the cross-variance only appears once the vessel is significantly deformed (the cross-variance is relatively flat around the plane of best focus for LSCI images), such that the two peaks in the kurtosis curve approximately coincide with the location where the cross-variance begins to appreciably increase. In addition, multiplication of kurtosis by cross-variance is approximately equivalent to removing the square in the denominator of Eq. (1), resulting in

3. Results

We performed two experiments to characterize the ability of ζ to better estimate the location of the plane of best focus compared to the SFI maps, showing that: (a) ζ is more sensitive to deviation from the plane of best focus and (b) ζ can be used as a more robust feedback measure than flow index in an autofocus routine.

3.1. Experimental setup

The experimental setup used in this work is shown in Fig. 3. The system consists of a microscope objective (Motic 10× Apo-plan, NA = 0.24) mounted on a motorized translation stage (Newport ESP300, Irvine, CA, USA). A f = 200 mm field lens is used to form an intermediate image which is subsequently projected onto a camera (QImaging Rolera EM-C², Surrey, BC, Canada) using a 4 − f imaging system. A spatial light modulator (SLM, Meadowlark Optics Inc., 512×512 Namatic Liquid Crystal, Frederick, CO, USA) is placed in the Fourier plane of the 4 − f system to serve as the active element for the autofocusing experiment described in section 3.3. The location of the plane of best focus is shifted by introducing a quadratic phase pattern into the beam with the SLM. This can be viewed as introducing a lens with variable lens power values at the SLM's surface location. A vertical-cavity surface-emitting laser (VCSEL, Vixar Inc., Plymouth, MN, USA) emitting coherent light at 680 nm was used in all experiments. A constant current source (Model 6220, Keithley Instruments, Cleveland, OH, USA) was used to drive the VCSEL with 3 mA constant current which resulted in a single mode of (high coherence) illumination.

 

Fig. 3 Experimental setup. An 10× objective lens and a f = 200 mm tube lens create an intermediate image, that is Fourier-transformed and projected onto an SLM. The SLM imposes quadratic phase patterns onto the beam, thus shifting the location of the plane of best focus. The image is then inverse-Fourier-transformed onto a camera. The SLM is designed to work only in one polarization, hence a polarizer is used to reject the polarization that would not be affected by the SLM.

Download Full Size | PDF

In all experiments we imaged the somatosensory cortex of an anaesthetized (2% isoflurane) healthy male Sprague Dawley rat (approx. 300 g) following a surgical protocol described in [20]. The animal study was performed in accordance with ethics protocols approved by the University of Toronto Animal Care Committee.

3.2. Kurtosis as a robust measure of misfocus

Figures 4(e-f) show the measured ζ and SFI values as a function of misfocus. The measurement is performed on a vessel portion selected from a field of view (Fig. 4(a)). To obtain this result, the microscope objective was translated vertically with respect to the vessel, across a span of 300 μm with 1 μm steps. The result shown indicates that for the particular vessel in question, ζ is both a less noisy (i.e. has lower frame-to-frame variance) measure of the position of the vessel with respect to best focus and provides a narrower localization of the location of best focus. We define the extent of localization of the range of best focus as a ‘range of ambiguity’, i.e. the range in z where there is uncertainty as to whether the system is aligned to the plane of best focus. The range of ambiguity of a quantity X(z) (either relative velocity or ζ), with a measurement standard deviation σ (X(z)), is defined as the span in z where $\overline{X(z)}+\sigma (X(z))>max[\overline{X(z)}-\sigma (X(z))]$. While our definition of the range of ambiguity in finding the plane of best focus is not unique, it provides a consistent quantitative means of comparison between different measures of misfocus by taking into account the measurement noise, σ(X). Combining the lower noise measurement with the narrower range of ambiguity of ζ, compared to SFI, allows us to obtain a 66 μm-wide window where ζ can indicate a change in misfocus while SFI can provide a 173 μm-wide window for the same vessel. The mean SFI values range between 167 [a.u.] to 148 [a.u.] for the range of ambiguity measured using ζ and between 167 [a.u.] to 130 [a.u.] for the range of ambiguity measured using LSCI. The FWHM of the flow profile, that is extracted from Fig. 4(d) and which is indicative of vessel diameter, ranges between 36.5 μm to 33.8 μm for the range of ambiguity measured using ζ, and between 40.1 μm to 32.1 μm for the range of ambiguity measured using LSCI. We stress here that the measurements shown in Fig. 4(e-f) have not been subjected to a temporal averaging filter and individual data points represent the distribution of single-frame measures. Figure 4(b-c) depicts the noise in single-frame SFI maps which can be reduced through temporal averaging. The result shown in Fig. 4(e-f) shows that the noise in a single-frame SFI map has far less impact on a measure of the SFI flow profile, ζ, than on the integrated SFI map as a whole. Consequently, using ζ enables a reduction in the noise associated with misfocus measurement without delaying the misfocus assessment response time as would be required using a temporal averaging approach to noise reduction.

 

Fig. 4 Comparison between the tranverse flow profile mean and the flow profile measure ζ as a means by which to estimate best focus location. (a) Field of view showing a SFI map of vessels in a somatosensory cortex in a rat, temporally averaged over 64 frames. The white square shows the ROI selected for analysis. (b) A zoomed-in depiction of the ROI selected from the temporally averaged SFI map in (a). (c) A non-averaged SFI map from a single speckle image in the selected ROI, showing the noise level of the SFI signal. (d) Transverse SFI flow profile vs. axial (z-) position of the objective with respect to the brain tissue. Every vertical slice represents a temporal average of 64 individually recentered SFI flow profiles acquired at a given depth. (e) Transverse flow profile mean (i.e., the average SFI from a recentered flow profile) verses z-position. The green curve shows the mean value of this measure for the 64 frames acquired at each z-step, while the shaded area represents the standard deviation. The vertical lines show the range of ambiguity in finding the best focus. (f) Our measure ζ as a function of misfocus (z-position). The vertical lines show the range of ambiguity in finding the best focus.

Download Full Size | PDF

3.3. Characterization of a closed-loop autofocus system for real-time correction of misfocus

A direct application of the measure ζ we propose is as a feedback mechanism in a LSCI auto-focus system. Robust, rapid autofocus for LSCI is challenging due to the inherent noise of the technique which manifests itself as large frame-to-frame variation in the flow index. Furthermore, the flow index values are sensitive to small-amplitude vibration which provide spurious flow information, e.g., as shown in reference [8].

We performed a controlled autofocus experiment to demonstrate the advantage of using our proposed measure compared to using the speckle flow index values directly. The objective lens was moved axially with respect to brain tissue, using the motorized translation stage, and the measure ζ was monitored over time, along with the SFI values. A simple threshold-based autofocus routine was implemented: the system was made to refocus using the SLM when ζ decreased below a threshold empirically set just below the average value of ζ at best focus. The refocusing routine maximized ζ using a simple binary search across a set of discrete focal length values, imposed through quadratic phase profiles on the SLM surface. In our particular implementation, we discretized the search space in z into 101 positions, which required the acquisition of 7 (= ceil[log₂ 101]) frames to find the location of ζ_max. The camera frame rate in our implementation was 10 fps, limited by the acquisition speed of the MATLAB Image Acquisition Toolbox. We note here that any number of focus search algorithms could have been used; we chose the binary search as a suitable compromise between simplicity of implementation and speed of focusing. While more sophisticated algorithms would yield faster convergence of the system to the plane of best focus, improving the speed of the focusing mechanism has not been one of our objectives in this study.

Figure 5 shows the results of our autofocus system for translations with amplitudes of 70 μm. Two advantages are inherently evident in using our measure. First, it shows significantly lower frame-to-frame noise (we did not use a moving average filter). As such, the low frame-to-frame noise allows us to robustly threshold the autofocus routine, which would not be possible using flow speed values directly obtained from Fig. 5(d). Secondly, since our proposed measure is based solely on the vessel shape rather than its flow values, it is not sensitive to changes in the flow values themselves (due to falling blood pressure or jerking movement in the system, for example). The effect of the correction of a single vessel on reducing the flow measurement error over the entire field of view is shown in Figs. 5(g-h). We define the measurement relative error as $|\overline{SFI(z)}-\overline{SFI(z=0)}|/\overline{SFI(z=0)}$, where $\overline{SFI(z)}$ is the flow index at a particular plane z, averaged over 100 frames, where z = 0 is the plane of best focus. The results in Fig. 5(g) show significant errors in flow index measurements, particularly along vessels’ edges, where uncorrected 70 μm axial movements (second and fourth tiles from the left) caused up to 75% error in the extracted flow speeds, observed through a significant change in color between adjacent tiles. Figure 5(h) shows that correcting for the axial displacement of a single vessel using ζ virtually eliminated these errors, as there is significantly less change in color between the four tiles in Fig. 5(h) compared to those in Fig. 5(g).

 

Fig. 5 Closed loop autofocus results using ζ as a feedback mechanism. (a) Full field of view relative flow velocity map of a rat cortex recorded on the camera. The white square shows the ROI selected for autofocus computation. (b) The axial position of the focal plane of the microscope objective with respect to the location of best focus. The shaded regions in (b-f) indicate that an attempt by the autofocus system to detect and correct for an induced misfocus. The unshaded regions indicate no attempt to correct for misfocus (i.e., flat phase mask applied to the SLM). (c) The value of the proposed metric ζ. The horizontal line shows the threshold for refocusing. (d) Speckle flow index inside the vessel. (e) A 15-frame moving average of the speckle flow index. (f) The power of a lens displayed on the SLM surface. (g) Axial motion-induced relative error across the entire field of view for axial objective positions of (from left to right) 0 μm, 70 μm, 0 μm, and −70 μm with respect to the plane of best focus. Every tile is an average of each focal condition for the uncorrected temporal windows. (h) Relative error, as in (g), with correction using the SLM.

Download Full Size | PDF

4. Discussion

Finding and maintaining the location of best focus in LSCI during an in-vivo imaging session is imperative for extracting quantitative results on flow speeds and vessel diameters. In this work we present a new image processing routine that quantifies optical defocusing using a statistical measure, ζ, and demonstrate its ability to better correct for misfocus in an autofocus system. Our proposed method’s key benefit is its more precise measurement of deviations from the plane of best focus compared to using the integrated SFI values directly. As demonstrated in Fig. 4, it provides a ×2.6 narrower localization of the vessel best focus compared to using the SFI directly which, in turn, results in ×2 more accurate reading of the SFI and ×2.9 more accurate reading of the flow profile width. Moreover, this measure is not sensitive to rapid motion artefacts. While such motion may cause potentially significant changes in the extracted flow speed values, the value of ζ would not show appreciable change since the vessel’s flow profile shape is not appreciably altered by uniform background movement.

Misfocus-related motion artefacts can also be reduced by lowering the numerical aperture of the system, which increases the systems depth of field. However, there are two key drawbacks to the above approach: 1) closing the aperture limits the light intensity incident on the camera which forces the user to work at higher exposure times, limiting the temporal resolution and 2) decreasing the numerical aperture increases the speckle size, limiting the spatial resolution.

We would like to note that we performed the autofocus routine on an amplitude of movement much greater than typical physiological movements due to breathing and pulsation in a rat cranial window. Indeed, for an appreciable change in SFI to be recorded a vessel’s misfocus needs to be comparable to the depth of field of the optical system. Hence, the amplitude of motion applied in Fig. 5(b) was chosen to illustrate our technique’s sensitivity rather than correct for a genuine defocusing artefact. We would like to further emphasize that although we used only a small-ROI in the field of view for the assessment of misfocus, we were able to correct for errors in the entire field of view (Fig. 5(g-h)). This is indicative of the robustness of our technique as only a single suitable vessel is required. Moreover, it emphasizes the utility of this approach for improving LSCI while only processing a small fraction of the acquired data.

The performance of our algorithm in accurately estimating misfocus depends upon choosing a vessel of an appropriate diameter. For example, a vessel whose diameter is larger than, or comparable to, the depth of field of the optical system will not deform significantly due to a small amount of misfocus and, therefore, the ability of the algorithm to correct for small deviation around the location of best focus would be limited. Conversely, a vessel that only spans several pixels in diameter does not allow for sufficient resolution in resolving the transverse flow profile and, as such, the algorithm would not produce a robust metric for misfocus. In our experience, choosing a vessel whose diameter is slightly smaller than the depth of field of the imaging system produces the best results.

Our proposed technique could be utilized in applications that require accurate measurements of location of best focus and/or vessel diameter. These applications include experiments aimed at measuring vasodilation and vasoconstriction, where changes of the order of <10% are typically measured [21]. Based on our results, for example, deviations from best focus of >52 μm would have >10% change in the flow profile FWHM. This implies that the range of ambiguity of ζ is sufficient for rapid tracking of vasodilation or vasoconstriction, while the range of ambiguity of the integrated SFI is not. Our technique could also benefit multi-modal imaging studies which combine Intrinsic Optical Signal Imaging (IOSI) and LSCI, where the LSCI modality can be used to ensure that the imaging plane remains in best focus for both modalities. Maintaining best focus for IOSI is particularly important since the changes in oxygenation that are extracted using IOSI are dependent solely on changes in contrast in the image (i.e. greylevel changes recorded on the camera) [3], and the degradation of contrast with misfocus would result in errors in extraction of quantitative oxygenation changes.

5. Conclusion

In this work we presented an image processing system and technique that enables sensitive monitoring of tissue location with respect to the plane of best focus by quantifying changes in vessel morphology away from best focus. Our proposed measure for misfocus, ζ, the product of flow profile kurtosis and the “apparent width” of the vessel, allows for 2.6-fold improvement in the sensitivity to misfocus as compared with the integrated speckle flow index maps, resulting in 2-fold better estimate of the relative flow speed. We demonstrated this measure as a feedback mechanism for an autofocus routine in a controlled experiment with large deviations from best focus, showing an ability to correct for movement with amplitude of 70 μm, using single-frame measurements in a fast autofocus routine. Future applications of this proposed method include accurate measurements of vessel diameters in experiments that study vasodilation and vasoconstriction, in multi-modal optical imaging, and in motion correction for in-vivo studies in awake animals.