## Abstract

Single-frame blood flow maps from laser speckle contrast imaging (LSCI) contain high spatiotemporal variation that obscures high spatial-frequency vascular features, making precise image registration for signal amplification challenging. In this work, novel bivariate standardized moment filters (BSMFs) were used to provide stable measures of vessel edge location, permitting a more robust LSCI registration. Relatedly, BSMFs enabled the stable reconstruction of vessel edges from sparsely distributed blood flow map outliers, which were found to retain most of the temporal dynamics. Consequently, data discarding and BSMF-based reconstruction enable efficient real-time quantitative LSCI data compression. Smaller LSCI-kernels produced log-normal blood flow distributions, enhancing sparse-to-dense inference.

© 2018 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 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. 4*f* correlator design: two 30 cm focal length lenses, ∅8 mm front and back apertures, optical wavelength 680 nm, Fresnel number 40, and 2048^{2} grid elements with 8 μm pitch. The Fresnel impulse response propogator of the complex field *U* over free-space distance *z* was

*x*is the spacing of the simulation grid (

*x*,

*y*) [24].

#### 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*〉

^{2}/[〈

*I*

^{2}〉 − 〈

*I*〉

^{2}]. Within the context of the two applications demonstrated here, the most useful post-

*k*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,

_{n}*ξ*, 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*. 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

_{xy}*ξ*

_{3}and

*ξ*

_{4}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

*ξ*

_{4}and

*ξ*

_{6}filters appear to provide better background feature suppression relative to edge-based features than the

*ξ*

_{3}and

*ξ*

_{5}filters. The

*ξ*

_{6}filter provided only slightly better feature visibility than the

*ξ*

_{4}filter. The similarity between these two filters implies that the

*ξ*

_{4}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)).

We quantified the SNR, $\u3008{\sigma}_{\text{signal}}^{2}\u3009/\u3008{\sigma}_{\text{noise}}^{2}\u3009={\u3008{\sigma}_{xy}^{2}\u3009}_{t}/{\u3008{\sigma}_{t}^{2}\u3009}_{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 *ξ*_{5} 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 *ξ*_{6} increasing the least. For this region, the *ξ*_{3} filter had a higher SNR than the *ξ*_{5} filter.

We found that BSMFs are able to enhance the accuracy of image displacement estimates. In particular, the *ξ*_{3} and *ξ*_{5} 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

*ξ*

_{4}

*previous frame*reference. The

*ξ*

_{3}and

*ξ*

_{5}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.

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 *ξ*_{3} and *ξ*_{5} 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 *ξ*_{4} and *ξ*_{6} filters also showed greater breathing related displacement specificity (see Fig. 2(b)). However, the *ξ*_{4} filter, in particular, showed the least specificity and was prone to vertical location estimation error.

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

*ξ*

_{4}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)).

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 *ξ*_{3} or *ξ*_{4} 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 *ξ*_{3} and *ξ*_{4} 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 *ξ*_{4} filter-directed registration reduced the accumulation of error, while the *ξ*_{3} 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 *ξ*_{3} and *ξ*_{4} 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).

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 *ξ*_{3} or *ξ*_{4} 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 *ξ*_{3} or *ξ*_{4} 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 *ξ*_{4} 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

*ξ*

_{4}filter map shown in Fig. 4(a) corresponds to square filter window with vertical value projection. This is the simplest and fastest

*ξ*

_{4}filter. As most vessel edges in this example are vertical this simple filter detects the majority of features. The

*ξ*

_{4}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 *ξ*_{4} (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 *ξ*_{4} (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 *ξ*_{4} (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*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,

_{n}*k*

_{3}, we see a log-normal distribution of blood flow values (excess kurtosis of −0.1 ± 2 × 10

^{−3}). For increasing

*k*, additional features begin to become evident as the distribution deviates from log-normal (excess kurtosis of −0.5, −0.7, & −0.8 ± 2 × 10

_{n}^{−3}for

*k*

_{5},

*k*

_{7}, &

*k*

_{9}, 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

^{−7}for

*k*

_{3}and +0.42 ± 2 × 10

^{−7}for

*k*

_{9}) and negative skew shift (−0.15 ± 1 × 10

^{−3}for

*k*

_{3}and −0.28 ± 1 × 10

^{−3}for

*k*

_{9}) 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*, yet the blood flow value distributions for small

_{n}*k*have less shape variability. The combined simplicity and shape stability of

_{n}*k*

_{3}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*distributions have greater similarity, with the most distinct distribution (

_{n}*k*

_{3}) 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*

^{2}> 0.99 ± 4 × 10

^{−5}, 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*

_{3}and

*k*

_{5}, 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*

_{3}and

*k*

_{5}range, respectively, in linear flow value space.

#### 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

_{γ}*ξ*

_{3}(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*map. Consequently,

_{n}*ξ*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*

_{5}kernel resulted in better arteriole signal retention, whereas the smaller

*k*

_{3}kernel resulted in more balanced signal retention. Note: bounds in correlation plots computed through transformation to the z-distributed parameter

*R′*= tanh

^{−1}

*R*, bounding ${R}_{\partial}^{\prime}={R}^{\prime}\pm {\sigma}_{{R}^{\prime}}$, and inversion ${R}_{\partial}={R}_{\partial}^{\prime \u2020}=\text{tanh}{R}_{\partial}^{\prime}$. The maximum

*ξ*

_{3}spatial frequency occurred at a higher level of sparsity for

*k*

_{3}than

*k*

_{5}but

*k*

_{3}had less spatial frequency content for pixel retention above 1 % (see Fig 6(e)).

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

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 *ξ*_{3} edge detection features, derived from *ξ*_{3} 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 *ξ*_{3} 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 *ξ*_{3} application on either sparse or dense single-frame blood flow maps. To quantify the stability and bias of edge features from *ξ*_{3} filters, vessel diameter estimates from *ξ*_{3} 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 *ξ*_{3} diameter estimate was based on the peak-to-peak distance from *ξ*_{3} 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/σ^{2} is rapidly lost. Conversely, the *ξ*_{3} based diameter estimates retain significant precision up to and including 1.0 ±.5 % pixel retention. Furthermore, for *ξ*_{3} both the pre-ictal and ictal vessel diameter estimates (*d*_{0} and *d*_{1}, 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 *ξ*_{3} estimate for assessing relative diameter change, $\frac{\mathrm{\Delta}d}{{d}_{0}}$, due to proportional *d*_{0} and *d*_{1} 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

*ξ*

_{3}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

*ξ*

_{3}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

*ξ*

_{3}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

*ξ*

_{3}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.

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*^{2} < 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.

## 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.

## 6. Appendix

## 6.1. Rolling sum algorithm for arbitrary directional moment calculation

For a given univariant data set, the most efficient calculation of a higher-order central moment is its defining formula:

*x*− 〈

*x*〉)

^{2}〉 = 〈

*x*

^{2}〉 − 〈

*x*〉

^{2}. The advantage of this formula is that if a new data point is to be included or an old one removed, the three summations over

*w*,

_{i}*w*·

_{i}*x*and ${w}_{i}\cdot {x}_{i}^{2}$ can be modified with a single addition/subtraction applied to each term, avoiding recalculation of the summations. We simply extend this approach to the bivariate central moments up to order

_{i}*γ*≤ 6; in terms of the centralized variable

*u*=

*x*− 〈

*x*〉 and

*v*=

*y*− 〈

*y*〉 we seek all moments, 〈

*u*〉, where 0 ≤

^{α}v^{β}*α*,

*β*and

*α*+

*β*≤

*γ*. These central moments can be computed for any group of pixels {

*x*,

*y*} from the summation terms,

*n*(

*n*+ 1)/2 terms. Section 6.2 presents the computational routine for finding all corresponding 〈

*u*〉 terms efficiently. The non-centralized sum terms Σ

^{α}v^{β}*can be updated with terms ${T}_{\alpha \beta}={w}_{i}\cdot {x}_{i}^{\alpha}\cdot {y}_{i}^{\beta}$ through addition and subtraction, Σ*

_{αβ}*= Σ*

_{αβ}*±*

_{αβ}*T*, permitting use within an efficient rolling sum algorithm. The standard rolling sum algorithm for a square window is shown in Fig. 10(a). The efficiency of the square window algorithm is independent of window size. The

_{αβ}*T*terms can be computed efficiently together using prior powers,

_{αβ}*T*=

_{αβ}*T*

_{α(β−1)}·

*y*, or,

_{i}*T*=

_{αβ}*T*

_{(α−1)β}·

*x*, where

_{i}*T*

_{00}=

*w*.

_{i}The 〈*u ^{α}v^{β}*〉 permit us to calculate the univariate higher moments in any orientation. Any new orientation can be defined by a linear combination of coordinates,

*w*=

*au*+

*bv*, where

*a*

^{2}+

*b*

^{2}= 1. The univariate moments 〈

*w*〉 with 1 ≤

^{γ}*γ*are calculated by taking the inner product of the “bivariant central moment vector” [〈

*u*〉, 〈

^{γ}*u*

^{γ−1}

*v*〉, . . . 〈

*v*

^{γ−1}

*u*〉, 〈

*v*〉 with an array of binomial expansion terms from (

^{γ}*a*+

*b*)

*as follows:*

^{γ}A circular sliding window was also employed to reduce the effect of window shape on moment orientation. The rolling sum algorithm we developed for a circular window is shown in Fig. 10(b). The efficiency of the circular algorithm is proportional to the window size as leading and lagging arcs are used for each window step. However, this algorithm is better suited to parallelization on modern processors than the square algorithm. A diagram of an exhaustive search of the five 4* ^{th}*-order moments assessed across the four possible

*π*/4 orientations for a circular window is shown in Fig. 10(c). Note that only odd sized windows (pixels across) were used to permit the superposition of resultant smaller filtered images over flow images during analysis.

The computational complexity of the brute force, square rolling sum, and circular rolling sum algorithms are *O*(*N*^{2}*b*^{2}), *O*(*N*^{2}), and *O*(*N*^{2}*b*), respectively, where *N* is the number of pixels in the image width/height and *b* is the width/height of the sliding window. The actual number of vector sum update operations per image is less than *N*^{2}*b*^{2}, 4*N*^{2}, and 2*N*^{2}*b*, respectively, for the three aforementioned algorithms. For a *ξ*_{4} filter with *π*/8 resolution and 51 pixel diameter applied to a 500 × 500 pixel image (such as in Fig. 4(c)), the rolling cirular rolling sum was > 800 × faster than using the non-computational standard moment formula with explicit pixel direction projection. Furthermore, the error introduced by catastrophic cancellation associated with the computational formulas was negligible, ${\sigma}_{\text{signal}}^{2}/{\sigma}_{\text{error}}^{2}=16\times {10}^{12}$. For sparse images, the regularization contribution to the circular rolling sum algorithm was pre-computed (*i.e.* all arc additions and subtractions). This enabled the majority of pixels to be accounted for with an efficiency exceeding the square rolling sum algorithm (we term this sparse acceleration).

## 6.2. Bivariant central moment calculation

The derivation of 〈*u ^{α}v^{β}*〉 terms from Σ

*terms defined by Eq. 2 is presented here. The 〈*

_{αβ}*u*〉 terms are used to calculate the directional moments 〈

^{α}v^{β}*w*〉 via Eq. 3. In this notation, the means are calculated from the rolling sums as follows $\u3008x\u3009={\mathrm{\Sigma}}_{10}\cdot {\mathrm{\Sigma}}_{00}^{-1}$ and $\u3008y\u3009={\mathrm{\Sigma}}_{01}\cdot {\mathrm{\Sigma}}_{00}^{-1}$. We then proceed to calculate the higher order standard moments 〈

^{γ}*u*〉 using lower order moments and the terms used in their computation. This includes the iteratively computed powers of the means, 〈

^{α}v^{β}*x*〉

*= 〈*

^{n}*x*〉

^{n−1}〈

*x*〉. Starting with the principal coordinate variance and covariance terms (2

*order),*

^{nd}*order terms which are the numerators in principal coordinate skew and co-skew.*

^{rd}*c*

_{21}= Σ

_{21}− 〈

*x*〉 · Σ

_{11}and

*c*

_{12}= Σ

_{12}− 〈

*y*〉 · Σ

_{11}. Similarly, we proceed to the 4

*order terms which are the numerators in principal coordinate kurtosis and co-kurtosis.*

^{th}*c*

_{31}= Σ

_{31}− 3 〈

*x*〉 ·

*c*

_{21}and

*c*

_{13}= Σ

_{13}− 3 〈

*y*〉 ·

*c*

_{12}. The 5th and 6th order terms can be calculated similarly.

## Funding

National Sciences and Engineering Research Council (NSERC) (RGPIN-355623-08); Canadian Institutes of Health Research (CIHR) (CPG-121050 and MOP-119603); Connaught Fund (498042); MITACS (IT02018).

## Acknowledgment

The authors would like to thank Iliya Sigal, Melanie A. Jeffery, Suzie Dufour and Hart Levy for their significant contributions to experimental studies from which data was derived for this proof-of-concept work.

## Disclosures

The authors declare that there are no conflicts of interest related to this article.

## References

**1. **D. Ringuette, M. A. Jeffrey, S. Dufour, P. L. Carlen, and O. Levi, “Continuous multi-modality brain imaging reveals modified neurovascular seizure response after intervention,” Biomed. Opt. Express **8**(2), 873–889 (2017). [CrossRef] [PubMed]

**2. **L. M. Richards, E. L. Towle, D. J. Fox, and A. K. Dunn, “Intraoperative laser speckle contrast imaging with retrospective motion correction for quantitative assessment of cerebral blood flow,” Neurophotonics **1**(1), 015006 (2014). [CrossRef]

**3. **S. E. Skipetrov, J. Peuser, R. Cerbino, P. Zakharov, B. Weber, and F. Scheffold, “Noise in laser speckle correlation and imaging techniques,” Opt. Express **18**(14), 14519–14534 (2010). [CrossRef] [PubMed]

**4. **I. Sigal, M. M Koletar, D. Ringuette, R. Gad, M. A. Jeffrey, P. L. Carlen, B. Stefanovic, and O. Levi, “Imaging brain activity during seizures in freely behaving rats using a miniature multi-modal imaging system,” Biomed. Opt. Express **7**(9), 3596–3609 (2016). [CrossRef] [PubMed]

**5. **R. Farraro, O. Fathi, and B. Choi, “Handheld, point-of-care laser speckle imaging,” J. Biomed. Opt. **21**(9), 094001 (2016). [CrossRef]

**6. **P. Miao, A. Rege, N. Li, N. V. Thakor, and S. Tong, “High resolution cerebral blood flow imaging by registered laser speckle contrast analysis,” IEEE Trans. Biomed. Eng. **57**(5), 1152–1157 (2010). [CrossRef] [PubMed]

**7. **S. J. Kirkpatrick, D. D. Duncan, and E. M. Wells-Gray, “Detrimental effects of speckle-pixel size matching in laser speckle contrast imaging,” Opt. Lett. **33**(24), 2886–2888 (2008). [CrossRef] [PubMed]

**8. **D. Coleman, “TIFF compression options: ZIP vs LZW” (TIFF compression options, 2018). https://havecamerawilltravel.com/photographer/tiff-image-compression/.

**9. **S. S Parikh, D. Ruiz, H. Kalva, G. Fernández-Escribano, and V. Adzic, “High bit-depth medical image compression with hevc,” IEEE J. Biomed. Health Informat. **22**(2), 552–560 (2018). [CrossRef]

**10. **K. Amolins, Y. Zhang, and P. Dare, “Wavelet based image fusion techniques-An introduction, review and comparison,” ISPRS J. Photogramm. Remote Sens. **62**(4), 249–263 (2007). [CrossRef]

**11. **D. Ringuette, I. Sigal, R. Gad, and O. Levi, “Reducing misfocus-related motion artefacts in laser speckle contrast imaging,” Biomed. Opt. Express **6**(1), 266–276 (2015). [CrossRef] [PubMed]

**12. **C.-H. Teh and R. T. Chin, “On image analysis by the methods of moments,” IEEE Trans. Pattern Anal. Mach. Intell. **10**(4), 496–513 (1988). [CrossRef]

**13. **D. Shen and H. H. S. Ip, “Discriminative wavelet shape descriptors for recognition of 2-d patterns,” Pattern Recognition **32**(2), 151–165 (1999). [CrossRef]

**14. **H. Levy, D. Ringuette, and O. Levi, “Rapid monitoring of cerebral ischemia dynamics using laser-based optical imaging of blood oxygenation and flow,” Biomed. Opt. Express , **3**(4), 777–791 (2012). [CrossRef] [PubMed]

**15. **P. Thevenaz, U. E Ruttimann, and M. Unser, “A pyramid approach to subpixel registration based on intensity,” IEEE Trans. Image Process. , **7**(1), 27–41 (1998). [CrossRef]

**16. **J. W. Goodman, “Speckle with a finite number of steps,” Appl. Opt. , **47**(4), A111–A118 (2008). [CrossRef] [PubMed]

**17. **I. D. Schizas and G. B. Giannakis, “Covariance eigenvector sparsity for compression and denoising,” IEEE Trans. Signal Process. **60**(5), 2408–2421 (2012). [CrossRef]

**18. **J. D. Shutler, M. S. Nixon, and C. J. Harris, “Global statistical description of temporal features,” Int. Arch. Photogramm. Remote Sens. **33**(B5/2; PART 5), 720–726 (2000).

**19. **J. D. Shutler and M. S. Nixon, “Zernike velocity moments for sequence-based description of moving features,” Image Vis. Comput. **24**(4), 343–356 (2006). [CrossRef]

**20. **D. Wen, Y. Jiang, H. Hua, R. Yu, Q. Gao, and Y. Zhang, “Laser speckle reduction based on compressive sensing and edge detection,” Proc. SPIE **8905**, 890506 (2013). [CrossRef]

**21. **S. M. Shams Kazmi, E. Faraji, M. A. Davis, Y. Huang, X. J. Zhang, and A. K. Dunn, “Flux or speed? examining speckle contrast imaging of vascular flows,” Biomed. Opt. Express **6**(7), 2588–2608 (2015). [CrossRef]

**22. **M. Hu, “Visual pattern recognition by moment invariants,” IRE Trans. Inf. Theory **8**(2), 179–187 (1962). [CrossRef]

**23. **J. F. Boyce and W. J. Hossack, “Moment invariants for pattern recognition,” Pattern Recognit. Lett. **1**(5–6), 451–456 (1983). [CrossRef]

**24. **D. G. Voelz, *Computational fourier optics: a MATLAB tutorial* (SPIE Press, 2011), Chap. 5.