1. Introduction

A technique for full-field blood flow imaging, laser speckle contrast imaging (LSCI), provides rapid snapshots of blood flow values. When averaged spatially over a suitably large region-of-interest (ROI), such as over a large vessel, relative blood flow values can be shown to have a high degree of accuracy and temporal precision [1]. Furthermore, when explicitly correcting for rapid temporal variations due to cardiac pulsation, the effective signal-to-noise ratio (SNR) can be significantly improved [1, 2]. Unfortunately, pseudo-random high-frequency spatial and temporal measurement variation burdens LSCI with a low SNR at small spatial and short temporal scales (as with all laser speckle correlation techniques [3]). To compensate, temporal averaging of spatially stable or accurately registered sequential images can be used to enhance the spatial resolution. However, the aforementioned variation in the LSCI blood flow maps (or unprocessed speckle images) is highly problematic for accurate sequential image registration as high contrast features are obscured. Consequently, potential short time-scale image alignment artifacts caused by animal breathing, specimen manipulation, or behavioral locomotion [4] are difficult to correct. Under ideal experimental conditions, relative movement can be limited to less than one micron. However, movement of tens-of-microns due to the aforementioned causes is not uncommon, particularly in cases where implant or restraint interfaces have become compromised. Consequently, LSCI image registration is often desired but faces unique challenges. Moreover, given the potential demonstrated for application of LSCI in many scenarios where movement is likely, such as intra-operatively [2] or point-of-care applications [5], there is significant benefit to robust registration for signal amplification. It has been shown that the perceived vessel visibility associated with macroscopic temporal speckle contrast maps are improved by prior speckle image registration based on persistent pixel-scale structural features [6]. However, it has not been demonstrated that robust LSCI registration is possible for microscopy-based LSCI, cortex-only field-of-views (FOVs), or large movements. Microscopy-based LSCI permits accuracy optimized (pixel-wise speckle contrast Nyquist criteria [7]) micro-vascular flowmetry with high optical collection efficiency and minimum speckle size relative to vessel size (higher resolution). The increased accuracy and resolution are associated with higher spatiotemporal signal variation and, as such, persistent pixel-scale structural features are absent. Moreover, for small FOVs, as extra-cortical features are absent, image registration must be directed solely by features extracted from tissue undergoing vascular dynamics. Furthermore, pixel-scale speckle-driven features are assuredly changed for large displacements.

In a pragmatically separate but fundamentally related issue, the quantitative high bit-depth images required by LSCI and other scientific imaging applications are difficult to compress. The on-disk size of such images cannot be reduced, and is often increased [8], by conventional compression schemes. Recent advancements and a characterization of the problem complexity associated with high-bit depth medical image compression are provided by reference [9]. This is a particularly difficult issue for LSCI where the high spatial frequency information prevents efficient compression with wavelet-based codecs, as wavelet transforms are effectively high pass filters [10]. Moreover, active image compression requires computational resources which can interfere with stable data transfer. Thus, due to the inherent noise of individual raw speckle images or single-frame blood flow maps, LSCI is limited both with respect to data compression efficiency and registration stability. In this work, an image filtering approach is introduced to mitigate these noise induced limitations.

We have previously demonstrated a sub-depth of field (DOF) active misfocus correction scheme for LSCI based on the fourth standardized moment of vessel cross-sections [11]. A covariance matrix corresponding to a circular ROI enabled vessel orientation estimation independent of lateral translation. Then translation insensitive profiles were computed by cropped rotation and projection of an inscribed square ROI with subsequent recursive one-dimensional cropping around the profile mean. Our misfocus measure, which we termed ζ, used a mixture of statistics from circular and square ROIs to provide a larger dynamic range for misfocus correction. It is well established that higher-order moments retain detailed image information but are sensitive to white noise [12,13], making them unsuitable for many image discrimination tasks [13]. Thus, we asserted that the sub-DOF behavior of ζ is due to vessel geometric features being closely associated with misfocus. However, within single-frame LSCI-derived blood flow maps, noise is multiplicative (ie., non-white noise). Consequently, the accuracy of ζ may be further due to spatial moment-based statistics in general having a high SNR when applied to LSCI-derived blood flow maps. Therefore, we are motivated to expand the characterization of spatial moment statistics and parameters affecting blood flow value distributions. Moreover, we seek to determine if spatial moment features can be further exploited for lateral motion correction in addition to axial. A practical limitation was the absence of an image filtering algorithm for ζ-like statistics.

In this work, we introduce an efficient rolling-sum approach for calculating central moment-based statistics across arbitrary spatial directions within a circular or square sliding window. We found that when applied to LSCI blood flow maps, bivariate standardized moment filters (BSMFs) produced sharp stable features associated with vessel edges. We demonstrate that these edge features can enhance lateral registration of sequential LSCI frames, thereby increasing the spatial resolution and accuracy of blood flow measurements. Furthermore, both temporal blood flow dynamics and spatial vessel edge dynamics could be reproduced from sparsely distributed outlying blood flow values, thereby permitting quantitative real-time data compression.

2. Methods

This work is a proof-of-concept demonstrating stable edge detection in LSCI blood flow maps. Stable edge-detection was exploited for image registration and data compression. We focused on imaging blood flow dynamics associated with disease models of seizure and stroke in superficial rodent cortex.

2.1. Imaging systems, surgical procedures, and speckle imaging simulation

Laser speckle images were collected on our previously published imaging systems to quantify the properties of BSMFs in LSCI. Three imaging systems were used: table-top [14], miniature head-mounted [4] and, primarily, a microscope integrated system [1]. The field-of-view of the systems was >2 mm, ∼2 mm and 200–400 μm, respectively. All three systems were optimized for concurrent LSCI and intrinsic optical signal imaging (IOSI) using coherence modulation of vertical-cavity surface-emitting lasers (VCSELs). The IOSI channel permits concurrent measurement of blood oxygenation spectroscopically. The LSCI illumination wavelength was 680 nm. The low spatial contrast and high temporal variation of low coherence near-infrared reflectance images make IOSI ill-suited to image registration (see Visualization 1).

The surgical procedure for the microscope integrated system [1] consisted of a craniotomy followed by the mounting of a perfusion system, which permitted objective lens immersion, illumination light-guide immersion, and topical convulsive drug (4-aminopyridine) application. The surgical procedure for the table-top imaging system [14] consisted of a craniotomy followed by an agarose stabilized cover-slip mount. Stroke was induced through carotid artery restriction via surgical suture. The miniature imaging system [4] required the implantation of a chronic optical window and device mount weeks before imaging. Absence seizure-like events were induced through intraperitoneal injection of a convulsive drug (pentylenetetrazol). The photothrombosis example (rose bengal 30 mg/kg i.v. with 545 ± 10 nm excitation) was conducted on the microscope integrated system without the immersion optics using a CD-1 mouse.

Fresnel propagation simulations of a microscope were used to assess location accuracy with known object perturbations. We simulated a 1 × mag. 4f correlator design: two 30 cm focal length lenses, ∅8 mm front and back apertures, optical wavelength 680 nm, Fresnel number 40, and 2048² grid elements with 8 μm pitch. The Fresnel impulse response propogator of the complex field U over free-space distance z was

2.2. Choice of edge enhancement filter

We investigated the properties of spatial moment-derived image filters applied to single-frame LSCI-based blood flow maps acquired using the brain imaging systems referenced in Section 2.1. Individual camera exposures of speckle intensity were converted to blood flow maps by applying an n × n pixel rolling sum speckle flow index filter, k_n = 〈I〉²/[〈I²〉 − 〈I〉²]. Within the context of the two applications demonstrated here, the most useful post-k_n image filters were those providing stable vessel edge detection. The parameter we found most useful was the maximum univariate standardized moment occurring in a set of directionally projected cross-sections. The projection was most efficiently calculated across arbitrary directions using bivariate centralized moments (see Appendix 6.1) which were efficiently calculated from the bivariate non-centralized moments (see Appendix 6.2) generated via a rolling-sum algorithm. For simplicity, we chose the term bivariate standardized moment filters (BSMFs) for all statistics derived in this manner. We introduce the symbol, ξ_γ, for the γ-order directionally optimized BSMF. Efficient BSMF algorithms were developed for both circular (with sparse acceleration) and square filter windows.

2.3. Filter directed lateral image registration

To improve the resolution and accuracy of LSCI blood flow maps, single-frame images were registered (translation and rotation) using the publicly available ImageJ plugin StackReg prior to temporal averaging. This plug-in provides a propagating registration where the previous frame anchors registration of the next. The algorithm is optimized for precise detection of sub-pixel movement between frames [15]. We trivially modified the algorithm such that sequential single-frame blood flow maps could be registered based on the registration computed for filtered images generated from these maps.

2.4. Sparsity enforcement scheme for data compression

The BSMFs were used for edge detection using fractions of the original data. We implemented a naive sparsity enforcement scheme in which only the top portion of blood flow values were retained and all values below a threshold were discarded. An encoding strategy was adopted such that the ratio of retained data (relative to an unprocessed frame) was twice the ratio of retained pixels: The retained 32-bit floating point flow values (threshold to max) were rescaled into unsigned 16-bit integers (0 to 65535). The retained pixel indexes were encoded as the row major distance between adjacent pixels (also unsigned 16-bit integers). For moment calculations, all blood flow values had a regularization offset of one.

3. Results

The results are presented as follows: 1) Low-noise edge-selective properties of BSMFs are characterized including the spatial precision associated with BSMFs for image displacement assessment. 2) Reduction of movement induced measurement error through sub-pixel registration is demonstrated, with quantification of registration precision enhancement. 3) The effect of window geometry, directional optimization and blood flow value rescaling on BSMFs are investigated. 4) The dependence of blood flow value distributions on flow index kernel choice is investigated. 5) The effect of explicit sparsity enforcement on fast temporal dynamics and BSMF-based edge reconstruction is presented. 6) The dynamic tracking of reconstructed edge features from sparse blood flow values during large flow perturbations for the accurate tracking of absolute vessel diameter is shown. 7) Lastly, the reconstruction of temporal dynamics associated with large flow perturbation is demonstrated through non-linear histogram re-normalization.

3.1. Low-noise edge-selective properties of bivariate standardized moment filters

We found a distinguishing feature of BSMFs, as applied to LSCI blood flow maps, is that they select precise features around vessel edges while retaining little noise from the original blood flow map. A computed blood flow map is inherently noisy and spatial averaging (or band-passing) retains a significant portion of this noise (see Fig. 1(a, e)). In fact, speckle suppression through both averaging and band-passing were found to increase registration error (data not shown). As we found previously [11], the orientation of an observed vessel can be stably computed from a covariance matrix, C_xy, corresponding to a circular window spanning a vessel. The longitudinal and cross-sectional directions correspond to eigenvectors of C_xy. The corresponding eigenvalues are extrema of variance ${\sigma }_{\text{max}}^2$ and ${\sigma }_{\text{min}}^2$, respectively. For smaller window diameters, the parameters ${\sigma }_{\text{max}}^2$ or ${\sigma }_{\text{min}}^2$ behave as edge selective filters (see Fig. 1(b, f)). However, the edge features are not highly discernible from background. Conversely, the BSMFs-derived parameters ξ₃ and ξ₄ were found to produce sharply defined features at vessel edges from single-frames (see Fig. 1(c,g)). The low specificity of the variance-based parameter in the denominator of all BSMFs did not appear to limit filter edge specificity. A similar phenomena was previously observed with respect to our sub-depth of field misfocus parameter ζ [11]. From observation of cross-sectional profiles, the ξ₄ and ξ₆ filters appear to provide better background feature suppression relative to edge-based features than the ξ₃ and ξ₅ filters. The ξ₆ filter provided only slightly better feature visibility than the ξ₄ filter. The similarity between these two filters implies that the ξ₄ filter is sufficient in most cases. Changing the size of the filter window results in a trade-off between background noise and the sharpness of edge features (see Fig. 1(d,h)).

 

Fig. 1: Stable edge detection properties of BSMFs. (a) Single-frame LSCI blood flow map. (b) Filters based on extrema of 2^nd-order moments i) ${\sigma }_{\text{min}}^2$ and ii) ${\sigma }_{\text{max}}^2$ applied to a. (c) The BSMFs i) ξ₃ and ii) ξ₄ applied to a. (d) ξ₄ as in c for i) smaller and ii) larger filter windows. Filter window sizes in b–d indicated by black circle. (e–h) Cross-section profiles from a–d, respectively, for dashed-line indicated in a. The red curve in e corresponds to vertical pixel averaging within white box in a. Higher-order BSMFs ξ₅ and ξ₆ also shown in g.

Download Full Size | PDF

We quantified the SNR, $〈{\sigma }_{\text{signal}}^2〉/〈{\sigma }_{\text{noise}}^2〉={〈{\sigma }_{xy}^2〉}_t/{〈{\sigma }_t^2〉}_{xy}$, of filtered images from the example shown above as spatial feature variation over temporal signal noise occurring during short (0.8 s) stable epochs between animal breaths (see Table 1). All BSMFs significantly improved the SNR relative to the original blood flow maps and second order filters. For the full FOV, the ξ₅ filter had the highest SNR and its temporally averaged filter (not shown) was the most specific to vessel edges. For a 100 ×100 μm region spanning the arteriole shown (for which the filter and flow values were high) the SNR was higher for all filters with ξ₆ increasing the least. For this region, the ξ₃ filter had a higher SNR than the ξ₅ filter.

Table 1:. The SNR associated with LSCI blood flow maps and associated BSMFs assessed during static epochs, for the full FOV in Fig. 1 and a 100 × 100 μm arteriole-centered ROI.

View Table | View all tables in this article

We found that BSMFs are able to enhance the accuracy of image displacement estimates. In particular, the ξ₃ and ξ₅ filters had a ∼ 2.3 to 2.9 fold higher intrinsic spatial precision than blood flow maps. We characterized the spatial stability of maps through assessment of mean location shifts, $r=\sqrt{x^2+y^2}$, implied by frame-to-frame cross-correlation during short (0.8 s) stable epochs between animal breaths (see Table 2). The full-frame cross-correlations were evaluated over ± 17 pixels (±13.6 μm in object plane). Two reference frame conditions were assessed: the first frame in an epoch and the previously occurring frame. For both conditions, the ${\sigma }_{\text{min}}^2$ and ${\sigma }_{\text{max}}^2$ filters decreased measured spatial stability. Conversely, the ξ_γ filters increased measured spatial stability, with the exception of the ξ₄ previous frame reference. The ξ₃ and ξ₅ filters had the least movement instability suggesting they are the filters best suited for the evaluation of image displacement. The average frame-to-frame error was lower for the previous frame registration condition for blood flow maps and all corresponding filters. This suggests there is an advantage to rolling or pyramidal registration approaches for LSCI signals.

Table 2:. The spatial location precision of LSCI blood flow maps and associated BSMFs assessed with cross-correlation during static epochs. Units are pixels (scale 0.8 μm/pixel).

View Table | View all tables in this article

Similarly, we used a cross-correlation to assess whether image displacements could be more accurately tracked by ξ filtered images. For the example chosen, there was a recurrent breathing related rightward displacement of 12 μm (15 pixels). The size and direction of the recurrent displacement was estimated visually a priori. The ξ₃ and ξ₅ filters predicted sample displacements consistent with the visual assessment of both artifact timing and magnitude (see Fig. 2(a)). Moreover, they showed higher specificity than blood flow maps which predicted similar displacements. The displacement estimates based on blood flow above the 17 pixel search range reflect error in vertical location estimation. The ξ₄ and ξ₆ filters also showed greater breathing related displacement specificity (see Fig. 2(b)). However, the ξ₄ filter, in particular, showed the least specificity and was prone to vertical location estimation error.

 

Fig. 2: Spatial precision assessment of filters via cross-correlation based estimation. (a) Comparison of inferred displacement from ξ₃ and ξ₅ filters against inference using blood flow maps. Actual displacement image is approximately 15 pixels. (b) Comparison of ξ₄ and ξ₆ filters for the same series in a. (c) Comparison of all ξ_γ for simulated imaging of an optical phase varying cross flow shape displaced 10 pixels at frame 20. Computed flow map and associated ξ₃ filter shown above graph. (d) Same as c for optical phase varying square.

Download Full Size | PDF

We performed Fresnel propagation simulations of an LSCI system as this allowed us to assess accuracy using known displacements. Furthermore, it provided a confirmation that our observed noise insensitivity is not solely dependent on some specific feature of our experimental procedure or system. Two different samples were assessed: 1) A cross of varying optical phase to reflect vessels with different orientations (see Fig. 2 (c)) and 2) a square of varying optical phase to provide a case without continuous edges (see Fig. 2 (d)). In both cases, the simulated phase varying object was displaced 10 pixels to the right. The relative accuracy of ξ_γ filters was consistent with the experimental results. In particular, the ξ₄ filter had the lowest displacement specificity and was distinctly worse for the more challenging square flow shape case.

3.2. Registration-based reduction of movement induced measurement errors

The stability associated with BSMFs was found to extend to a more precise sub-pixel registration approach [15]. The chosen algorithm used nonlinear least-square intensity difference minimization with pyramid (coarse-to-fine) spline fitting. This algorithm was the one we found in practice to be most reliable for registering short sequences of LSCI blood flow maps. However, its accuracy was reduced by large sample motions and gradual accumulation of alignment error. Here we show that sub-pixel registration error was reduced using ξ_γ filter-direct registration.

In the example from Section 3.1, high resolution blood flow maps can be generated from the LSCI images occurring between the breathing related motion artifacts (see Fig. 3(a)). Conversely, blood flow maps acquired during lateral sample motion have distorted vessel edges (see Fig. 3(b)).

 

Fig. 3: Improved sub-pixel registration of blood flow maps through BSMF guidance. (a) A high resolution blood flow map produced by temporally averaging maps occurring between breathing related motions (generated from fifteen 0.8 s stable epochs within a 16 s epoch). All graphs b–d, f, & g share a common correspondence with respect to registration strategy: stable reference (blue), ξ₃ filter directed (red), ξ₄ filter-directed (green), no registration (orange), and self-directed (black). (b) Vessel-edge cross-sections for motion artifact associated blood flow maps. The paired (same color) curves represent the cross-section distributions (μ ± σ) for fourteen blood flow maps each associated with a 0.3 s breath-related motion within the same 16 s epoch. (c) Artery adjacent region blood flow time trace distributions associated with the same motion artifacts in b (each registration was initiated from the indicated time zero). (d) Correlation loss of single-frame blood flow maps with reference map a. Distributions correspond to repeat series registrations each initialed at 1 of 5 sequential motion artifacts. (e) Correction for motion artifact due to a 34 s unintentional lateral stage drift. (f) Cross-sections corresponding to region indicated in e. (g) Correction for 11 s lateral motion artifact due to stroke induction. (h) Cross-sections for central vessel in g.

Download Full Size | PDF

 

Fig. 4: Effect of filter window shape, orientation searching and non-linear value re-scale. (a) Map of square window ξ₄ filter based on only vertical univariate projections. Extension of filter in a to include (b) eight univariate orientations, and (c) a circular filter window. (d & e) Logarithmic value rescaling prior application of ξ₄ in b & c, respectively. (e) Same as c with prior logarithmic rescaling. (f) Same as c with prior exponential rescaling.

Download Full Size | PDF

For these short motion artifacts, high resolution temporally averaged blood flow maps can be generated from such image sequences after the application of the sub-pixel registration algorithm directly to single-frame flow maps (i.e. self-directed registration). Using either ξ₃ or ξ₄ filter-directed registration increased the accuracy of high resolution flow profile estimations. This can be seen as the narrower distribution width, μ ± σ, of location specific flow measurements at the vessel boundary (see Table 3(b)). The higher accuracy of ξ₃ and ξ₄ filter-directed registration is most clearly observable from the consistency of blood flow time traces before and after correction of this recurrent motion artifact (see Fig. 3(c)). Conversely, self-directed sub-pixel registration cannot consistently remove motion artifacts, observable as a 3.6 ± .3× increase in signal variation post artifact (see Table 3(c)). Moreover, errors in self-directed registration accumulate more rapidly over time (see Fig. 3(d)). This is prohibitive for sequential frame registration approaches (such as this algorithm) which benefit from coherence driven similarity in temporally adjacent frames and can account for large feature changes, provided they occur gradually. Using ξ₄ filter-directed registration reduced the accumulation of error, while the ξ₃ filter almost completely suppressed the accumulation of error (see Table 3(d)). In the presense of changing edge features during vessel dilation (i.e. onset of seizure-like event shown in Fig. 7&9(a)), the ξ₃ and ξ₄ filter-directed registration stabilities converged and self-directed registration became only slightly more stable (see Table 3(d′)). Conversely, exaggerating the rate of vessel diameter change by synthetically interleaving more distinct maps could eliminate the filter-based stability advantage (data not shown). As such, fast LSCI acquistion is likely necessary for stable non-elastic image registration (i.e. for image alignment using only rotation and translation).

Table 3:. Accuracy of registration approaches. Corresponding to Fig. 3: (b) The variation of edge profiles during movement artifact. (c) The variation of time traces before and after registration of recurrent artifact. (d) Rate of correlation loss with respect to static reference.

View Table | View all tables in this article

During one imaging experiment, there was a temporary lateral subject displacement due to unintentional physical contact with the microscope translation stage. For this larger (∼ 45 μm) and longer (∼ 34 s) motion artifact, self-directed registration was worse than no registration at all (see Fig. 3(e, f)). Conversely, using ξ₃ or ξ₄ filter-directed registration produced a corrected imaging series for which the temporal average was highly correlated with a blood flow map generated from frames acquired prior to movement. For a similar motion artifact (∼ 11 s and ∼ 50 μm) over a larger FOV (2 mm), self registration performance was superior to no registration (see Fig. 3(g, h)). However, using ξ₃ or ξ₄ filter-directed registration again produced a temporally averaged blood flow map that more closely matched an associated static reference map. This small motion was due to suture tension applied to the carotid artery for induction of ischemia.

We found that ξ_γ filters within about 20% of the estimated optimum filter size would produce comparable registration. In general, attempts to suppress noise such as the application of spatial low-pass filters or a rolling average only made registration less stable leading to rapid walk-off (data not shown). This is likely due to a relative absence of distinct features in the image.

3.3. Effect of filter window shape, direction searching, flow value re-scaling

In this section we look at the interaction effects of perturbation from optimal BSMF parameters. In addition to providing an intuition for accuracy/efficiency trade-offs, features leading to stabilization of filters are revealed. Three modifications are investigated: 1) Filter window shape; square filter windows are more efficient but projections along diagonal orientations are different than horizontal/vertical. 2) Projection orientation searching; more comprehensive direction searching would be expected to provide more precise features. 3) Blood flow map value re-scaling; suppression of outlying values should, in general, enhance edge detection. We chose the ξ₄ filter for this characterization as its lower SNR and sharp features should be more susceptible to perturbation.

In the proceeding section, we focus solely on the properties of BSMFs based on circular filter windows applied directly to blood flow maps. The computational efficiency of ξ_γ filters depends on filter window shape and resolution of direction searching (see Appendix 6.1). The ξ₄ filter map shown in Fig. 4(a) corresponds to square filter window with vertical value projection. This is the simplest and fastest ξ₄ filter. As most vessel edges in this example are vertical this simple filter detects the majority of features. The ξ₄ based on eight orientations, as shown in Fig. 4(b), highlighted additional edges but introduced non-specific background features. Conversely, using a circular window, as shown in Fig. 4(c), highlights even further vessel edges without the background signal from the parenchyma. Consequently, outside of this characterization, section circular windows were used exclusively. For the blood flow maps we investigated, using optimally selected filter diameters, the ξ_γ filters sharpness was negligibly improved beyond π/8 projection angle resolution.

To investigate the contribution of outlying flow values, we perturbed the measured blood flow values with a non-linear rescaling (either logarithmic or exponential). A logarithmic rescaling was applied before applying a square window ξ₄ (see Fig. 4(d)). This significantly reduced edge feature continuity and increased background signal. A logarithmic rescaling was applied before applying the circular window prior to ξ₄ (see Fig. 4(e)). This also increased background noise but to a lesser extent. These results suggest that circular windows are more robust to perturbation. Moreover, reducing the influence of outliers reduces filter feature specificity. As outlying flow values are most perturbed by logarithmic rescaling, this suggests outliers directly contribute to feature stabilization. The converse re-scaling (i.e. exponential) was applied before applying the circular window ξ₄ (see Fig. 4(f)). Individual filter maps were more varied frame-to-frame but median filter behavior shown was consistent with vessel edges and had low non-specific parenchymal signal. This suggested that outlying flow values support BSMF-based edge reconstruction but redistributing these and intermediate values increases signal variation. Consequently, we investigated the influence of outlying flow values more directly in Section 3.5.

3.4. Smaller speckle flow index kernels produce simpler distributions of sampled blood flow values

We found the sparsity of a single-frame blood flow map is dependent on the size of the flow index kernel, k_n. Smaller k_n have a greater range of outliers and simpler distribution shapes (most easily seen in a logarithmic scale, see Fig. 5(a)). In this case, for the limit of the smallest odd kernel, k₃, we see a log-normal distribution of blood flow values (excess kurtosis of −0.1 ± 2 × 10⁻³). For increasing k_n, additional features begin to become evident as the distribution deviates from log-normal (excess kurtosis of −0.5, −0.7, & −0.8 ± 2 × 10⁻³ for k₅, k₇, & k₉, respectively). These features include a second flow maximum appearing below the central flow maximum and a subtle bulging along the upper value tail. Increased blood flow due to a seizure-like event results in a positive mean shift (+0.43 ± 2 × 10⁻⁷ for k₃ and +0.42 ± 2 × 10⁻⁷ for k₉) and negative skew shift (−0.15 ± 1 × 10⁻³ for k₃ and −0.28 ± 1 × 10⁻³ for k₉) in the flow value histogram (see Fig. 5(b)). Based on the aforementioned measures, calculated changes in relative blood flow are similar for large and small k_n, yet the blood flow value distributions for small k_n have less shape variability. The combined simplicity and shape stability of k₃ implies its distribution can be most easily inferred given a subset of distribution samples. The lack of complexity would also imply that information is lost. However, after temporal averaging, the flow value distributions are far less dependent on kernel size (see Fig. 5(c)). Specifically, all four k_n distributions have greater similarity, with the most distinct distribution (k₃) appearing to have a constant offset (i.e. scaling factor not effecting relative flow calculations) and less sharp local extrema. Moreover, the predicted relative distribution shifts are still similar across kernel size (see Fig. 5(d)) and the temporally averaged frames from different kernel sizes are all highly correlated (R² > 0.99 ± 4 × 10⁻⁵, assuming spatial low-passing of the largest kernel diameter to remove resolution differences). Consequently, smaller kernel sizes appear to benefit from distribution simplicity without losing information sought in lower-noise temporally averaged maps. Therefore, we choose to focus on the smaller kernels, k₃ and k₅, for our expanded investigation on the influence of sparsity enforcement. Note: in Fig. 5(a) the upper 1% of values in flow distribution corresponds to 97 % and 78 % of the k₃ and k₅ range, respectively, in linear flow value space.

 

Fig. 5: Distribution of blood flow index values based on kernel size k_n. (a) Distribution for first four odd kernels n = 3, 5, 7, &9. (b) Distribution shift associated with ictal event for n = 3&9. (c & d) Distributions after temporal averaging of frames in a & b, respectively.

Download Full Size | PDF

3.5. Retention of the edge feature detection properties of BSMFs and subtle rapid temporal dynamics in the presence of sparsity enforcement

The vessel edge selection properties of ξ_γ filters were preserved when applied to images in which only sampled flow values above a threshold were retained. Sampled flow values above a sufficient threshold were pseudo-sparsely distributed (highest values often adjacent) within a single-frame and their locations were highly variable frame-to-frame. The spatial features observed in high resolution blood flow maps (generated by temporal averaging maps within short static 0.8 s epochs) degenerate rapidly with increased sparsity (see Fig. 6(a, e)). Conversely, the disappearance of edge-related features is more gradual for single-frame ξ_γ filtered images, such as ξ₃ (see Fig. 6(b, e)). In fact, edge features first become sharper with increasing sparsity, eventually approximating those derived from dense maps at a given sparsity. Moreover, as the edge features disappear with further sparsity enforcement, the filtered maps approximate the shape of the vessel dense k_n map. Consequently, ξ_γ filters enable spatial feature reconstruction from a small fraction of the sampled flow values. Furthermore, the simplicity of flow value distributions and their transformations, particularly for smaller kernel sizes (see Section 3.4), enables the highest flow values to retain the fast temporal dynamics associated with cardiac pulsation (see Fig. 6(c,d)). In this example, a consistent arteriole (lead) and venous (lag) phase relationship was observable for 0.25 % or 10^−2.6 pixel retention. Choosing the larger k₅ kernel resulted in better arteriole signal retention, whereas the smaller k₃ kernel resulted in more balanced signal retention. Note: bounds in correlation plots computed through transformation to the z-distributed parameter R′ = tanh⁻¹ R, bounding $R_{\partial }^{\prime }=R^{\prime }\pm {\sigma }_{R^{\prime }}$, and inversion $R_{\partial }=R_{\partial }^{\prime †}=\text{tanh}{R}_{\partial }^{\prime }$. The maximum ξ₃ spatial frequency occurred at a higher level of sparsity for k₃ than k₅ but k₃ had less spatial frequency content for pixel retention above 1 % (see Fig 6(e)).

 

Fig. 6: Spatial reconstruction and fast temporal dynamic preservation in the presence of sparsity enforcement. (a) Static epoch blood flow maps (25 frame average) with reduced pixel retention over five orders of magnitude (factors of 10^−1/5). (b) The ξ₃ filter applied to single-frame “sparse” blood flow maps from a (25 frame average). All maps in a and b were spatially low-pass filtered at the kernel size, k₅, for visual clarity. (c) Temporal flow profile of artery and vein for 25 frame static epoch during four cardiac cycle events (occurring at 5 Hz) for i) 10^−1.6, ii) 10^−2.6, iii) 10^−3.6 and iv) 10^−4.6, pixel retention (dashed curves represent no sparsity condition). (d) Correlation, R², of sparse and dense temporal profiles from c across all retention levels shown in a for both i) the k₅ results in a and ii) for k₃ (dashed curves are [(R′ ± σ_R′)^†]²). (e) Sum of spatial Fourier coefficients (features > 6 pixels) for both k₃ and k₅ and their corresponding ξ₃ filters (all curves are μ ± σ).

Download Full Size | PDF

We have shown that these filters enable spatiotemporal feature reconstruction from a trivial sparsity enforcement scheme. These filters close the loop between spatial and temporal retention of underlying vascular features in the presence of sparsity. However, it must be noted that sparsity enforcement reduces the stability of filters frame-to-frame. Consequently, registration and sparse reconstruction cannot be trivially combined. In the remaining results sections, we demonstrate that sparse encodings can be used to reconstruct vascular spatiotemporal dynamics associated with physiological events.

3.6. Edge detection properties of BSMF are robust to large flow and diameter changes enabling vessel shape and diameter tracking from sparse data

 

Fig. 7: Measuring vasodilation associated with seizure from sparse blood flow maps. (a–c) Superimposed pre-ictal (green) and ictal (red) maps of i) blood flow, ii) sparse blood flow, and iii) ξ₃ of single-capture sparse blood flow maps. All maps were produced from temporal average within single static epoch of 0.8–0.9 s. (d) Cross-sections for region in a (mean projection along the short axis of the indicated rectangle). (e) Cross-sections across sparsity for regions in a–c. All cross-section d–e regions match the ξ₃ kernel width to account for filter-based averaging. (f) Relative vessel diameter (pre-ictal no sparsity normalized to one) based on the estimated edge locations across sparsity from the cross sections in e. The edge features based on ξ₃ are contrasted with edge location estimation based on the FWHM. Paired same-color curves correspond to μ ± σ for 5 sequential static epochs. (g) The pre-ictal versus ictal difference in relative vessel diameters for example a including diameter estimation based on profile standard deviation.

Download Full Size | PDF

To determine if the retained structural information in sparsely sampled flow values is useful for detecting changes in vascular structure, we compared the sparse structural feature reconstruction before and during seizure-like event (ie, pre-ictal v. ictal). Such events involve changes in arteriole diameter and increased background capillary flow due to neurovascular coupling, both of which may bias edge localization. The artery flow profiles before dilation (green) and at peak dilation (red) are shown in Fig. 7(a–d)–i. The corresponding sparse blood flow maps, with 1 % of the sampled values retained, have poorly defined noisy edges for which precise edge estimation appears difficult (see Fig. 7(a–d)–ii). Conversely, the ξ₃ edge detection features, derived from ξ₃ application on single-frame sparse blood flow maps, appear to precisely mirror the vessel diameters before and during seizure-like events (see Fig. 7(a–d)–iii). Moreover, while the sparse flow profiles continually get narrower with increased sparsity, the ξ₃ profiles first widen then more gradually narrow and appear to match the dense filter profiles at 1.0 ±.5 % pixel retention (see Fig. 7(e)). Consequently, similar edge detection features can be derived from ξ₃ application on either sparse or dense single-frame blood flow maps. To quantify the stability and bias of edge features from ξ₃ filters, vessel diameter estimates from ξ₃ maps were contrasted with diameter estimates from their corresponding flow cross-section (see Fig. 7(f)). The FWHM was chosen for direct flow-map-based vessel diameter estimation as it is also based on edge approximation and provides accurate diameter measurements within dense blood flow maps without assuming a specific profile shape. The ξ₃ diameter estimate was based on the peak-to-peak distance from ξ₃ map cross-sections. We found the relative diameter change, $\frac{\mathrm{\Delta }d}{d_0}$, due to seizure was similar for both estimates from the dense maps. However, with increasing sparsity, the FWHM underestimates the absolute diameter, d, and the difference in diameter, Δd. Moreover, measurement precision, 1/σ² is rapidly lost. Conversely, the ξ₃ based diameter estimates retain significant precision up to and including 1.0 ±.5 % pixel retention. Furthermore, for ξ₃ both the pre-ictal and ictal vessel diameter estimates (d₀ and d₁, respectively) reverse an initial diameter overestimation, resulting in the 1.0 ±.5 % pixel retention derived estimates closely matching the dense derived estimates. We found that vessel diameter estimation based on the flow profile standard deviation underestimated absolute diameter with increased sparsity similar to the FWHM (data not shown). However, this measure was comparable to the ξ₃ estimate for assessing relative diameter change, $\frac{\mathrm{\Delta }d}{d_0}$, due to proportional d₀ and d₁ underestimation (see Fig. 7(g)). This indicates that sparse blood flow maps inherently contain the structural information, some of which can be assessed with second order moment-based statistics. Both estimates for the change in relative diameter were accurate for any pixel retention fractions at or above 0.40 ± .19 %. However, only ξ_γ filters were found to permit precise location and orientation estimation of vessel edge features from sparse blood flow maps for absolute diameter assessment. Furthermore, reconstructed edge features help guide the pre-processing stages of absolute diameter computation: 1) vessel region/segment selection, 2) region orienting for projection, 3) cross-section profile cropping, and 4) correction for vessel curvature. An example of the structural flexibility of cortical vessels is shown in Fig. 7(b)–iii where vessel curvature is higher when dilated. These precise location and orientation features also improve our ability to register slower displacement artifacts in sparse data as sparsity enforcement produces noise in unfiltered temporally averaged blood flow maps analogous to single-frame maps.

The arbitrary structural feature reconstruction achievable from sparse blood flow maps with ξ_γ filtering can be more easily seen by looking at examples of perturbations over larger FOVs. For instance, during photo-thrombotic ischemia, blood flow can become restricted within vessels (see Fig. 8(a)). In this case, a winding blood flow pattern is observed in the large vein in the center of the FOV. The stability of the ξ₃ filtering approach allows this pattern to be clearly reconstructed from sparse blood flow maps. To capture a large range of vessel sizes over the larger FOV, a small ξ₃ kernel was used. The rate of convergence to stable edge features in the presence of sparsity is proportional to kernel diameter. Consequently, additional frames beyond a static epoch were required to produce this map. A second example of feature reconstruction using a large FOV is provided by tracking of an absence-like seizure with our miniature device (see Fig. 8(b)). Using a ξ₃ filter optimized for large vessels, the dynamics of the full FOV are reproduced and a vasodilation can be observed in a large vessel. Enforcing sparsity on a smaller ROI and applying a smaller ξ₃ filter enabled visualization of small vessel dilation. A third example over a larger FOV is provided by imaging the application of a global ischemic stroke model with an incidental ischemia-induced spreading depression (see Fig. 8(c, d)). The decrease in vessel diameter post-stroke (due to restricted blood flow) and further decrease upon spreading depression (likely due to compensatory flow at the depression onset location) can be seen in the edge filter cross-sections. Lastly, an example of concurrent vasodilation and vasoconstriction is provided by the monitoring of a waking animal using our miniature imaging device (see Fig. 8(e)). This example emphasizes the need for reconstruction of individual vessel dynamics.

 

Fig. 8: Edge detection of features within a larger FOV. (a) Photothrombotic stroke induced through fluorescent excitation of rose bengal. i) Stroke (green) and pre-stroke (red) blood flow maps. ii) Corresponding temporally averaged 1% sparsity-enforced single-capture applied ξ₃ filter maps. (b) Vasodilation associated with absence-like seizure imaged with miniature imaging device. Panels i & ii same as a but for pre-ictal (green) and ictal (red). Panel iii is application of smaller ξ₃ for sparsity enforcement on region neglecting largest vessels. (c) Global stroke through carotid artery obstruction. i) Blood flow maps of pre-stroke (green), stroke (red), and ischemia-induced spreading depression (blue). ii) Corresponding temporally averaged 1% sparsity-enforced single-capture applied ξ₄ filter maps. (d) Cross section for middle vessel from c (yellow box). (e) Miniature device imaging of concurrent vessel dilation and constriction during animal waking; anesthetized (green) and awake (red).

Download Full Size | PDF

It should be noted that ξ_γ filters do not depend on sparse outliers (i.e. they are robust to subtle image smoothing) but are consistent despite them. Moreover, one can exploit their comparable behavior in the presence or absence of complete flow map data for structural feature reconstruction in the latter case.

3.7. Preservation of temporal dynamics associated with large flow perturbation in presence of sparsity

The vasodilation observed during seizure-like events is associated with a large increase in blood flow across the whole FOV. This is due to the hyper-synchronous activity inducing dilation around a large population of neuro-vascular units. The outlying blood flow values, in addition to enabling reconstruction of structural changes associated with seizure-like events, were found to enable reconstruction of the slower (≤ 1 Hz) time-scale temporal dynamics associated with such large flow increases (see Fig. 9(a)). The generation of accurate blood flow traces required a simple correction for thresholding-induced non-linear up-shifting of computed distribution means. The flow value histogram associated with baseline blood flow maps was scaled across increments of ± 5 % to derive an interpolated non-linear lookup table for the distribution mean inferred from the mean of values above a threshold. Without this non-linear correction, the onset and termination rates are underestimated with respect to peak flow and overestimated with respect to baseline. To use a signal variable one-to-one inference, pixel retention was defined from baseline and the associated threshold was applied throughout the series, as opposed to prior sections where pixel retention was fixed per static epoch. As such, for this scheme the pixel retention fraction is higher during periods of increased blood flow. Consequently, we also investigated a decreasing flow trend associated with global ischemic stroke (see Fig. 9(b)). For both the example of flow increase (due to seizure) and decrease (due to ischemia), the re-scaled sparse flow traces had low dynamical error (1 − R² < 3 %) and absolute error (normalized RMSE < 10 %) with respect to the dense blood flow traces down to 0.10 ±.05 % pixel retention (see Table 4, only major orders shown). As the histogram is a log-normal distribution, this interpolated relationship could, in principle, also have been estimated without ever acquiring dense baseline values. Additionally, a bivariate non-linear mean interpolation, accounting for a dynamic threshold, could permit a fixed pixel retention factor, as would be required for stable data transfer when system bandwidth is running near its limit.

 

Fig. 9: Temporal dynamics associated large blood flow changes are preserved in sparse values. The dependency of full-field blood flow temporal dynamics on sparsity during a (a) seizure-like event and (b) ischemic stroke. In both instances, level of sparsity set by threshold calibrated to baseline.

Download Full Size | PDF

Table 4:. The absolute error (normalized RMSE) and dynamical error (1 − R²) of sparse reconstructed temporal flow traces against dense flow trace. Note: the paired RMSE values are baseline and peak seizure/stroke flow increase/decrease in Fig 9, respectively.

View Table | View all tables in this article

4. Discussion

We demonstrated that BSMFs provide stable edge detection for noisy LSCI-derived single-capture blood flow maps (see Section 3.1). These BSMFs were initially investigated as we previously found the related ζ-statistics could provide sub-DOF misfocus localization [11]. Furthermore, BSMF could be modified into a rolling sum algorithm (see Appendix 6.1) giving them a computational advantage over other filter choices, such as Gabor filters (which require explicit convolution). We demonstrated that blood flow map-applied BSMF maps had a higher spatial feature SNR and location specificity than their source blood flow maps. We exploited this enhanced location specificity to improve the accuracy of sub-pixel registration applied to LSCI, thereby removing error in blood flow maps and time series (see Section 3.2). Enhanced registration was also demonstrated during vasodilation, when edge features were elastically changing, as the acquistion speed of LSCI enables rolling non-elastic alignment of similar feature maps. The characterization of different BSMF geometries and blood flow value re-scaling suggested that sparsely distributed outlying flow values contributed to, rather than detracted from, edge detection stability (see Section 3.3). We found that choosing smaller speckle flow index kernels increased flow value sparsity and decreased distribution complexity, without affecting relative flow calculations associated with physiological changes (see Section 3.4). Through sparsity enforcement on blood flow maps we found that from only the outliers of these distributions we could reconstruct vessel edge features and the temporal pulsation dynamics of individual vessels (Section 3.5). Moreover, these detected edge features were robust to large flow perturbations enabling the tracking of structural dynamics (see Section 3.6). Furthermore, the simplicity of the blood flow value distributions enabled accurate inference of the relative blood flow dynamics associated with these large changes using only outlying values (Section 3.7). Consequently, sparsity enforcement enables efficient LSCI data compression suitable for real-time implementation even with limited computational resources.

Speckle patterns have an approximately negative exponential intensity distribution [16], which is not particularly sparse. The observed sparsity in the blood flow maps is a property of the highly non-linear response of the speckle flow index calculation to image smoothness (functional analysis not shown). Here we demonstrated that the multiplicative noise associated with speckle flow index maps does not mask underlying structural features when assessed with BSMFs. We assert that the most likely source of edge feature stability, despite the high spatiotemporal variation of outliers, is that the bivariate mean reference point mitigates outlier variation within a given window. The unintuitive component of our result is that biasing toward outliers affects feature noise but not spatial convergence. Within the field of image processing, moment based denoising (with covariance eigenvector decomposition) has been demonstrated for non-biomedical images [17]. The results of Ref. [17] vary distinctly from those presented here, and bear little resemblance to the edge-based feature detection scheme outlined in this work. However, their results highlight a continuity to the literature for standard moments exploiting sparsity to account for noise. In other contexts, moments have been shown to be useful for velocimetry compression of moving object image sequences [18,19]. Others have applied the principles of compressive sensing-based edge detection for speckle reduction [20]. Their approach appears to produce a low-passed equivalent of speckle reduction. Conceptually, these two aforementioned works, Refs. [17] and [20], are the closest approaches we have found related to the structural feature extraction and sparse reconstruction scheme presented here, respectively.

A limitation of the registration approach demonstrated here is that it requires distinct flow differences, precluding application during severe ischemia, through skull (non-adaptive) imaging, or in cases of severe misfocus. There are also several limitations with respect to our data compression strategy: 1) the parenchymal blood flow values are under-represented, 2) closed-loop updating of the threshold is required to maintain the flow measurement dynamic range and/or stable read-out, and 3) for some specimens or implementations, flow value distribution may remain too complex for outlier guided inference. These issues can be partially mitigated through: 1) region adjusted thresholding, 2) concurrent retention of low-spatial resolution temporal dynamics, and/or 3) incremental threshold application based on application specific criteria. However, all of these solutions reduce data compression efficiency.

There are several imminently feasible extensions of this work: 1) Combining active multi-vessel ζ-driven active axial correction with post-hoc ξ-driven later correction and electrocardiography (ECG)-driven pulsation correction would exploit the full dynamic range achieveble with LSCI. This should promote studies of subtle neurovascular coupling (NVC) changes with reduced dependence on trial averaging for signal amplification. 2) Accurate segmentation of concurrent IOSI and LSCI series using BSMF-based edge features derived from individual LSCI captures. Spatially averaging within segmented regions should provide both signal amplification and efficient storage for both temporal blood oxygenation and flowmetry data. Moreover, segmentation parameters, such as spline curves, could be used for more efficient storage of spatial features. 3) Calibration of speckle flow index versus absolute velocity using accurate referencing of lateral motion artifacts. For any given LSCI measurement, there is an ambiguity between blood velocity and total flux [21] for which calibration could help determine the proportional contribution in any given region. 4) Adaptive optics permit coherent imaging using specimen proximate optics other than lenses (such as optical diffusers, multi-mode fiber optics and even superficial tissue) by correcting for the transmission matrix associated with arbitrary linear optics. Low-noise BSMF edge features could improve the per capture LSCI information content leading to more efficient phase optimization, thereby expanding the range of application under which LSCI can be implemented. 5) Using rotationally-invariant bivariate moment statistics [22, 23] to generate vessel branching pattern clusterings for automated ROI comparison or larger scale comparisons of vascular morphology. The circular rolling sum algorithm developed here permits the efficient calculation of these parameters which should be accurate for sparse single-frame captures enabling efficient vascular pattern storage.

5. Conclusion

We demonstrated BSMFs enable accurate LSCI spatial registration and provide a robust mechanism for vessel edge feature reconstruction from sparsely distributed outlying blood flow values. Third order BSMFs provided a 4.5 − 5.9× improvement in SNR and 2.6 − 2.9× improvement in location precision and reduced sub-pixel registration error accumulation by 23 ± 5×. The properties of BSMFs lead to the discovery that the highest blood flow values contain much of the spatiotemporal dynamics contained in LSCI. Accurate tracking of relative vessel diameter and temporal blood flow dynamics was demonstrated using only the top 0.40 ± .05 % and 0.10 ± .05 % of flow values, respectively. Accurate absolute vessel diameter estimation typically required retention of only the top 1.0 ± .5 % of blood flow values. Consequently, due to the low computational cost of flow calculations and threshold application, active sparsity enforcement and BSMF-based feature reconstruction enable a compressed sensing-like scheme for LSCI data acquisition.

 

Fig. 10: Rolling sum algorithms and direction searching. (a) Rolling sum algorithm for square sliding window. Sub-columns of window height are produced for all terms Σ_αβ in the first row. The first window is initialized by adding sub-columns and the window is moved horizontally by adding and removing sub-columns. The sub-columns are moved vertically by adding and subtracting pixels. Every pixel and sub-column term is added once and subtracted, at most, once. (b) Rolling sum algorithm for a circular sliding window. A circular window is first initialized and proceeds though successive arc additions and subtractions to trace the entire image. Every pixel is added and subtracted once times the circle diameter. (c) Exhaustive search across directions using centered expectation values.

Download Full Size | PDF