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

U(x,y;z)=FFT1{FFT{eikz(iλz)1eik2z(x2+y2)}(Δx)2FFT{U(x,y;0)}}
where FFT {·} is the 2D fast Fourier transform including zero-frequency centering and Δ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, kn = 〈I2/[〈I2〉 − 〈I2]. Within the context of the two applications demonstrated here, the most useful post-kn 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, Cxy, corresponding to a circular window spanning a vessel. The longitudinal and cross-sectional directions correspond to eigenvectors of Cxy. The corresponding eigenvalues are extrema of variance σmax2 and σmin2, respectively. For smaller window diameters, the parameters σmax2 or σmin2 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 ξ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)).

 

Fig. 1: Stable edge detection properties of BSMFs. (a) Single-frame LSCI blood flow map. (b) Filters based on extrema of 2nd-order moments i) σmin2 and ii) σmax2 applied to a. (c) The BSMFs i) ξ3 and ii) ξ4 applied to a. (d) ξ4 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 ξ5 and ξ6 also shown in g.

Download Full Size | PPT Slide | PDF

We quantified the SNR, σsignal2/σnoise2=σxy2t/σt2xy, 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.

Tables Icon

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.

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=x2+y2, 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 σmin2 and σmax2 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.

Tables Icon

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

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.

 

Fig. 2: Spatial precision assessment of filters via cross-correlation based estimation. (a) Comparison of inferred displacement from ξ3 and ξ5 filters against inference using blood flow maps. Actual displacement image is approximately 15 pixels. (b) Comparison of ξ4 and ξ6 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 ξ3 filter shown above graph. (d) Same as c for optical phase varying square.

Download Full Size | PPT Slide | 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 ξ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)).

 

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), ξ3 filter directed (red), ξ4 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 | PPT Slide | PDF

 

Fig. 4: Effect of filter window shape, orientation searching and non-linear value re-scale. (a) Map of square window ξ4 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 ξ4 in b & c, respectively. (e) Same as c with prior logarithmic rescaling. (f) Same as c with prior exponential rescaling.

Download Full Size | PPT Slide | 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 ξ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).

Tables Icon

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.

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, kn. Smaller kn 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, k3, we see a log-normal distribution of blood flow values (excess kurtosis of −0.1 ± 2 × 10−3). For increasing kn, additional features begin to become evident as the distribution deviates from log-normal (excess kurtosis of −0.5, −0.7, & −0.8 ± 2 × 10−3 for k5, k7, & k9, 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 k3 and +0.42 ± 2 × 10−7 for k9) and negative skew shift (−0.15 ± 1 × 10−3 for k3 and −0.28 ± 1 × 10−3 for k9) 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 kn, yet the blood flow value distributions for small kn have less shape variability. The combined simplicity and shape stability of k3 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 kn distributions have greater similarity, with the most distinct distribution (k3) 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 (R2 > 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, k3 and k5, 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 k3 and k5 range, respectively, in linear flow value space.

 

Fig. 5: Distribution of blood flow index values based on kernel size kn. (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 | PPT Slide | 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 ξ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 kn 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 k5 kernel resulted in better arteriole signal retention, whereas the smaller k3 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=R±σR, and inversion R=R=tanhR. The maximum ξ3 spatial frequency occurred at a higher level of sparsity for k3 than k5 but k3 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 ξ3 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, k5, 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, R2, of sparse and dense temporal profiles from c across all retention levels shown in a for both i) the k5 results in a and ii) for k3 (dashed curves are [(R′ ± σR′)]2). (e) Sum of spatial Fourier coefficients (features > 6 pixels) for both k3 and k5 and their corresponding ξ3 filters (all curves are μ ± σ).

Download Full Size | PPT Slide | 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) ξ3 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 ξ3 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 ξ3 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 | PPT Slide | 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 ξ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, Δdd0, 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 (d0 and d1, 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, Δdd0, due to proportional d0 and d1 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.

 

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 ξ3 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 ξ3 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 ξ4 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 | PPT Slide | 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 − R2 < 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 | PPT Slide | PDF

Tables Icon

Table 4:. The absolute error (normalized RMSE) and dynamical error (1 − R2) 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.

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 | PPT Slide | PDF

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:

(xxα)=i=0nwi(xix)αi=0nwi,wherex=i=0nwixii=0nwi
For the case of the second-order central moment (ie., variance) it can be calculated with the well known computational formula, 〈(x − 〈x〉)2〉 = 〈x2〉 − 〈x2. 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 wi, wi · xi and wixi2 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 γ ≤ 6; in terms of the centralized variable u = x − 〈x〉 and v = y − 〈y〉 we seek all moments, 〈uαvβ〉, where 0 ≤ α, β and α + βγ. These central moments can be computed for any group of pixels {x, y} from the summation terms,
Σαβ=i=0nwixiαyiβ
for which there are n (n + 1)/2 terms. Section 6.2 presents the computational routine for finding all corresponding 〈uαvβ〉 terms efficiently. The non-centralized sum terms Σαβ can be updated with terms Tαβ=wixiαyiβ 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) · yi, or, Tαβ = T(α−1)β · xi, where T00 = wi.

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 a2 + b2 = 1. The univariate moments 〈wγ〉 with 1 ≤ γ are calculated by taking the inner product of the “bivariant central moment vector” [〈uγ〉, 〈uγ−1v〉, . . . 〈vγ−1u〉, 〈vγ〉 with an array of binomial expansion terms from (a + b)γ as follows:

wγ=i=0γ(γi)aγibiuγivi
This directional projection equation utilizes the same principle exploited in the derivation of rotationally invariant moment-based measures [22].

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 4th-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(N2b2), O(N2), and O(N2b), 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 N2b2, 4N2, and 2N2b, 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, σsignal2/σerror2=16×1012. 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αvβ〉 terms are used to calculate the directional moments 〈wγ〉 via Eq. 3. In this notation, the means are calculated from the rolling sums as follows x=Σ10Σ001 and y=Σ01Σ001. We then proceed to calculate the higher order standard moments 〈uαvβ〉 using lower order moments and the terms used in their computation. This includes the iteratively computed powers of the means, 〈xn = 〈xn−1x〉. Starting with the principal coordinate variance and covariance terms (2nd order),

u2=Σ20Σ001x2,uv=Σ11Σ001xy,v2=Σ02Σ001y2
then proceeding to the 3rd order terms which are the numerators in principal coordinate skew and co-skew.
u3=Σ30Σ0013xu2x3,u2v=c21Σ001xuvyu2v3=Σ03Σ0013yv2y3,v2u=c21Σ001yuvxv2
where c21 = Σ21 − 〈x〉 · Σ11 and c12 = Σ12 − 〈y〉 · Σ11. Similarly, we proceed to the 4th order terms which are the numerators in principal coordinate kurtosis and co-kurtosis.
u4=Σ40Σ0014xu36x2u2x4,u3v=c31Σ001y(x3+u3)v4=Σ04Σ0014yv36y2v2y4,v3u=c13Σ001x(y3+v3)u2v2=(Σ222yc212xc12)Σ001+x2v2+y2(u2x2)
where c31 = Σ31 − 3 〈x〉 · c21 and c13 = Σ13 − 3 〈y〉 · c12. 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.

References

  • View by:
  • |
  • |
  • |

  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.

2018 (1)

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]

2017 (1)

2016 (2)

2015 (2)

2014 (1)

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]

2013 (1)

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]

2012 (2)

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]

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

2010 (2)

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]

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]

2008 (2)

2007 (1)

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]

2006 (1)

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]

2000 (1)

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

1999 (1)

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]

1998 (1)

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]

1988 (1)

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]

1983 (1)

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

1962 (1)

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

Adzic, V.

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]

Amolins, K.

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]

Boyce, J. F.

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

Carlen, P. L.

Cerbino, R.

Chin, R. T.

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]

Choi, B.

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

Dare, P.

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]

Davis, M. A.

Dufour, S.

Duncan, D. D.

Dunn, A. K.

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]

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]

Faraji, E.

Farraro, R.

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

Fathi, O.

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

Fernández-Escribano, G.

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]

Fox, D. J.

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]

Gad, R.

Gao, Q.

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]

Giannakis, G. B.

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

Goodman, J. W.

Harris, C. J.

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

Hossack, W. J.

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

Hu, M.

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

Hua, H.

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]

Huang, Y.

Ip, H. H. S.

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]

Jeffrey, M. A.

Jiang, Y.

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]

Kalva, H.

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]

Kirkpatrick, S. J.

Koletar, M. M

Levi, O.

Levy, H.

Li, N.

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]

Miao, P.

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]

Nixon, M. S.

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]

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

Parikh, S. S

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]

Peuser, J.

Rege, A.

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]

Richards, L. M.

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]

Ringuette, D.

Ruiz, D.

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]

Ruttimann, U. E

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]

Scheffold, F.

Schizas, I. D.

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

Shams Kazmi, S. M.

Shen, D.

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]

Shutler, J. D.

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]

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

Sigal, I.

Skipetrov, S. E.

Stefanovic, B.

Teh, C.-H.

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]

Thakor, N. V.

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]

Thevenaz, P.

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]

Tong, S.

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]

Towle, E. L.

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]

Unser, M.

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]

Voelz, D. G.

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

Weber, B.

Wells-Gray, E. M.

Wen, D.

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]

Yu, R.

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]

Zakharov, P.

Zhang, X. J.

Zhang, Y.

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]

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]

Appl. Opt. (1)

Biomed. Opt. Express (5)

IEEE J. Biomed. Health Informat. (1)

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]

IEEE Trans. Biomed. Eng. (1)

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]

IEEE Trans. Image Process. (1)

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]

IEEE Trans. Pattern Anal. Mach. Intell. (1)

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]

IEEE Trans. Signal Process. (1)

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

Image Vis. Comput. (1)

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]

Int. Arch. Photogramm. Remote Sens. (1)

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

IRE Trans. Inf. Theory (1)

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

ISPRS J. Photogramm. Remote Sens. (1)

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]

J. Biomed. Opt. (1)

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

Neurophotonics (1)

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]

Opt. Express (1)

Opt. Lett. (1)

Pattern Recognit. Lett. (1)

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

Pattern Recognition (1)

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]

Proc. SPIE (1)

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]

Other (2)

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

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

Supplementary Material (1)

NameDescription
» Visualization 1       Laser speckle flow index (SFI) maps and low coherence near-infrared (NIR) reflectance images. The SFI maps have distinct vascular features obscured by small-scale rapid variation. Conversely, NIR reflectance images have little tissue/vessel contras

Cited By

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

Alert me when this article is cited.


Figures (10)

Fig. 1:
Fig. 1: Stable edge detection properties of BSMFs. (a) Single-frame LSCI blood flow map. (b) Filters based on extrema of 2nd-order moments i) σ min 2 and ii) σ max 2 applied to a. (c) The BSMFs i) ξ3 and ii) ξ4 applied to a. (d) ξ4 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 ξ5 and ξ6 also shown in g.
Fig. 2:
Fig. 2: Spatial precision assessment of filters via cross-correlation based estimation. (a) Comparison of inferred displacement from ξ3 and ξ5 filters against inference using blood flow maps. Actual displacement image is approximately 15 pixels. (b) Comparison of ξ4 and ξ6 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 ξ3 filter shown above graph. (d) Same as c for optical phase varying square.
Fig. 3:
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), ξ3 filter directed (red), ξ4 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.
Fig. 4:
Fig. 4: Effect of filter window shape, orientation searching and non-linear value re-scale. (a) Map of square window ξ4 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 ξ4 in b & c, respectively. (e) Same as c with prior logarithmic rescaling. (f) Same as c with prior exponential rescaling.
Fig. 5:
Fig. 5: Distribution of blood flow index values based on kernel size kn. (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.
Fig. 6:
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 ξ3 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, k5, 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, R2, of sparse and dense temporal profiles from c across all retention levels shown in a for both i) the k5 results in a and ii) for k3 (dashed curves are [(R′ ± σR′)]2). (e) Sum of spatial Fourier coefficients (features > 6 pixels) for both k3 and k5 and their corresponding ξ3 filters (all curves are μ ± σ).
Fig. 7:
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) ξ3 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 ξ3 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 ξ3 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.
Fig. 8:
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 ξ3 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 ξ3 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 ξ4 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).
Fig. 9:
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.
Fig. 10:
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.

Tables (4)

Tables Icon

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.

Tables Icon

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

Tables Icon

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.

Tables Icon

Table 4: The absolute error (normalized RMSE) and dynamical error (1 − R2) 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.

Equations (7)

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

U ( x , y ; z ) = FFT 1 { FFT { e i k z ( i λ z ) 1 e i k 2 z ( x 2 + y 2 ) } ( Δ x ) 2 FFT { U ( x , y ; 0 ) } }
( x x α ) = i = 0 n w i ( x i x ) α i = 0 n w i , where x = i = 0 n w i x i i = 0 n w i
Σ α β = i = 0 n w i x i α y i β
w γ = i = 0 γ ( γ i ) a γ i b i u γ i v i
u 2 = Σ 20 Σ 00 1 x 2 , u v = Σ 11 Σ 00 1 x y , v 2 = Σ 02 Σ 00 1 y 2
u 3 = Σ 30 Σ 00 1 3 x u 2 x 3 , u 2 v = c 21 Σ 00 1 x u v y u 2 v 3 = Σ 03 Σ 00 1 3 y v 2 y 3 , v 2 u = c 21 Σ 00 1 y u v x v 2
u 4 = Σ 40 Σ 00 1 4 x u 3 6 x 2 u 2 x 4 , u 3 v = c 31 Σ 00 1 y ( x 3 + u 3 ) v 4 = Σ 04 Σ 00 1 4 y v 3 6 y 2 v 2 y 4 , v 3 u = c 13 Σ 00 1 x ( y 3 + v 3 ) u 2 v 2 = ( Σ 22 2 y c 21 2 x c 12 ) Σ 00 1 + x 2 v 2 + y 2 ( u 2 x 2 )

Metrics