Improved blood velocity measurements with a hybrid image filtering and iterative Radon transform algorithm

Neural activity leads to hemodynamic changes which can be detected by functional magnetic resonance imaging (fMRI). The determination of blood flow changes in individual vessels is an important aspect of understanding these hemodynamic signals. Blood flow can be calculated from the measurements of vessel diameter and blood velocity. When using line-scan imaging, the movement of blood in the vessel leads to streaks in space-time images, where streak angle is a function of the blood velocity. A variety of methods have been proposed to determine blood velocity from such space-time image sequences. Of these, the Radon transform is relatively easy to implement and has fast data processing. However, the precision of the velocity measurements is dependent on the number of Radon transforms performed, which creates a trade-off between the processing speed and measurement precision. In addition, factors like image contrast, imaging depth, image acquisition speed, and movement artifacts especially in large mammals, can potentially lead to data acquisition that results in erroneous velocity measurements. Here we show that pre-processing the data with a Sobel filter and iterative application of Radon transforms address these issues and provide more accurate blood velocity measurements. Improved signal quality of the image as a result of Sobel filtering increases the accuracy and the iterative Radon transform offers both increased precision and an order of magnitude faster implementation of velocity measurements. This algorithm does not use a priori knowledge of angle information and therefore is sensitive to sudden changes in blood flow. It can be applied on any set of space-time images with red blood cell (RBC) streaks, commonly acquired through line-scan imaging or reconstructed from full-frame, time-lapse images of the vasculature.


INTRODUCTION
The need to understand the relationship between neuronal and vascular activity has become crucial with the rise in the use of functional magnetic resonance imaging (fMRI) in research (Smith, 2012). Two-photon and other imaging techniques like bright field microscopy have provided us with the means to follow changes in neuronal and vascular activity at high spatial and temporal resolution. Neuronal activity can be determined by imaging calcium signals as a surrogate for spiking activity (Tsien, 1988;Stosiek et al., 2003;Ohki et al., 2005;Tian et al., 2009;Akerboom et al., 2013), while blood flow can be calculated from the vessel diameter and blood velocity.
Vessel diameter can be measured by labeling the vessel wall, for instance with an artery-specific dye , or by labeling the vessel lumen with circulating fluorescein-dextran (Schaffer et al., 2006;Drew et al., 2011). Blood velocity can be measured from the speed of red blood cells (RBCs) inside the vessel by labeling either the plasma or RBCs. Intravascular injection of large dextran-conjugated fluorescent dyes lead to labeling of the plasma. Alternatively, RBC-labeling with smaller lipophilic fluorescent dyes can also be used for studying vascular blood flow (Kamoun et al., 2010). These selective plasma or RBC-labeling methods lead to high contrast between RBCs and plasma, resulting in streaks of RBCs in the space-time images of vessels ( Figures 1A,B). A given streak angle depends on the distance traveled by the RBC per unit time. Therefore, blood velocity can be determined from the slope of such streaks.
Calculation of the slope of the RBC streaks has been proposed using many different algorithms including line fitting (Zhang et al., 2005), singular value decomposition (SVD) (Kleinfeld et al., 1998), Radon transform (Drew et al., 2010), and Fourier transform (Autio et al., 2011;Kim et al., 2012). Calculation of velocity using global energy minimization has also been proposed (Deneux et al., 2011). The Radon transform is one of the most powerful blood velocity measurement algorithms because of its easy implementation, accuracy, speed, and robustness. Radon transform calculates the slope from the angle that provides the maximum variance in the projections of the line-scan image at a specified angle space. However, it can still produce erroneous results when the line-scan images have low contrast and/or include artifacts (Figures 1C,D). Also, the precision of the velocity measurements is based on the angle step-size used for the Radon transform. Thus, selecting the number of angle steps to process with the Radon transform involves a trade-off between computation time and the precision of the measured streak angle.
Here we propose an enhanced algorithm that increases the accuracy and speed of velocity measurements: Sobel filtering of the image is followed by iterative application of multiple Radon transforms with increasing angle precision. Contrary to vertical (temporal) or whole-image demeaning, Sobel filtering (Sobel, 1978) efficiently removes the baseline or slow trends from the line-scan image. This translates into more accurate angle (and therefore velocity) measurements from line-scan images. Our development of an iterative Radon transform is built on the idea of adaptive application of Radon transforms with finer angle steps (see Drew et al., 2010 and our Materials and Methods). The first iteration of Radon transforms uses large angle steps across the full range of angles with progressively smaller steps in the subsequent iterations, allowing very precise angle calculations without the trade-off of exceedingly long computation times.

EXPERIMENTAL METHODS
Long Evans rats (postnatal days 29-48) and cats (postnatal day 28-adult) of both sexes were used to collect images from the primary visual cortex. The Institutional Animal Care and Use Committee of the Medical University of South Carolina approved all animal experiments. Detailed experimental setup, anesthesia, surgical procedures, and drifting grating visual stimulation protocols are described elsewhere Shen et al., 2012). Briefly, rats were anesthetized with intraperitoneal injection of a fentanyl cocktail (0.04-0.06 mg kg −1 fentanyl citrate, 3.75-6.25 mg kg −1 midazolam, and 0.19-0.31 mg kg −1 dexmedetomidine). Cats were anesthetized with 1-2% isoflurane. After making a 2-3 mm craniotomy and removal of dura, the vasculature was visualized with circulating 2000 kDa Fluorescein Dextran (FD2000S, Sigma), 75 kDa Texas Red Dextran (Life Technologies), or Alexa Fluor 633 Hydrazide (Life Technologies). Imaging was performed with a custom-designed two-photon microscope and a waterimmersion objective (Olympus XLUMPLFLN 20×, 1.0 NA or Olympus LUMPlanFI/IR 40×, 0.8 NA). Circulating RBCs were visualized as dark disks against the fluorescent-labeled plasma. These RBCs appear as dark bands or streaks in space-time images. In both rats and cats, movements of the brain were minimized by sealing the craniotomy with agarose and a coverglass O'Herron et al., 2012). In cats, a lumbar suspension and pneumothorax were needed to further reduce the brain movements of cardiac and respiratory origin. Typically, brain movement was limited to 1-2 μm in rats and kittens. But in some adult cats, movements of up to 5 μm were present.

DATA COLLECTION PARAMETERS
Fractional change in blood velocity is best detected when the RBC streak-angle is ∼45 • (see Equations 16,17). However, it is not always possible to acquire line-scan images with RBC streaks at this angle because of fast flow and/or relatively slow imaging, e.g., when using galvanometers. In such cases, the line-scan acquisition rate can be increased by imaging a smaller portion of the vessel, and/or by reducing the spatial resolution and pixel dwell-time of the scan, as long as the RBC-plasma contrast is sufficiently high. Additionally, the end-to-start flyback time during scanning can be reduced by combining a low-magnification objective with a compensatory zoom, e.g., 20× objective at zoom 4 instead of 40× objective at zoom 2. We observed that this slight advantage in speed gave superior line-scan acquisitions with a steeper slope of the RBC streaks and more reliable velocity measurements, e.g., in 30-50 μm diameter arterioles that have high blood speeds. Additionally, flyback speeds can be improved by driving the beam steering hardware to the fastest possible speeds.
For a vessel with blood velocity of v mm/s captured with linescan parameters of x μm/pixel spatial resolution and t ms/line temporal resolution (line-scan time), the angle of the streak θ is, The value of θ depends on blood velocity as well as spatial and temporal resolution of the line-scan image. If changing the imaging parameters offers a new angle θ new as a result of new spatial resolution of x new μm/pixel and temporal resolution of t new ms/line, the relationship between the new and old angles can be shown as, Where k is defined as, With the current streak angle θ and the desired streak angle θ new , the imaging speed factor k can be calculated as, In order to achieve the preferred angle of ∼45 • , the image acquisition speed should be changed by a factor of k = tan(θ). If the line-scan imaging parameters that give 60 • streak angle can be changed to improve the acquisition speeds by 30%, i.e., k = 1.3, the resultant streaks will have an angle of 53.11 • (Equation 2). A streak angle of 45 • can be achieved if the same imaging can be made faster by ∼73% (k = 1.73, Equation 4). The simplest way to change the image acquisition speeds (or t, line-scan time) is by changing the spatial resolution, pixel dwell time, and/or decreasing the pixel length of an imaged vessel segment. The caveat here is that the RBC streaks must maintain sufficient contrast to provide a signal-to-noise ratio (SNR) that enables further data processing and velocity measurements on the image (Figure 2).

SOBEL FILTERING OF THE LINE-SCAN IMAGES
Line-scan images often need to be filtered to remove artifacts that create luminance variations unrelated to the flow of RBCs. Artifacts may be static or slowly changing with time (relative to the line-scan duration) such as the artifacts caused by breathing or heartbeats in larger animals. Filtering is commonly performed by temporal demeaning of these space-time images, which is achieved by subtracting the mean value of all the pixels on the time axis at a given point in space. However, temporal demeaning can only address the time-invariant artifacts. Sobel filtering of the line-scan images efficiently removes both types of artifacts. We used Sobel filtering in the time-dimension of the linescan images (vertical Sobel filter, since our images have line-scans stacked vertically). The Sobel filter used here is 3 × 3 pixels and is applied by the convolution of the line-scan image I (x, t) with the vertical Sobel operator (S) as, Such filtering efficiently enhances horizontal edges and reduces vertical edges (Figure 3). This minimizes slow time-varying components of the line-scan images while sparing the RBC streak  edges. We did not find it necessary to use a Sobel operator with a larger kernel size, e.g., 5 × 5 or 7 × 7. Two-photon resolution produces sharp RBC streak edges, the thickness of which does not span more than a pixel or two. Moreover, line-to-line differences in the contrast of streak edges are best enhanced when the kernel height is 3.

Radon transform
The Radon transform of line-scan images has been described in detail earlier (Drew et al., 2010). Briefly, the Radon transform F (r, θ) of the image F (x, t) at a given angle θ can be defined as, F (τ sin(θ)+r cos(θ), −τ cos(θ)+r sin(θ)) dτ (7) The angle of the RBC streaks in the line-scan image is determined by the angle θ max at which the maximum variance is found, Blood velocity (v) is then calculated by scaling the tangent of this angle by the temporal ( t) and spatial ( One such Radon transform (angles spaced at 1 • ) is shown in Figure 4.

Effect of pixel resolution of the image on the precision of velocity detection
The Radon transform algorithm determines velocity from spacetime images by detecting an angle with maximum variance in the Radon transform space of the sampled angle values. While the Radon transform space can be sampled at any set of angles, the angle steps finer than the precision offered by the image only add to the computational load. The angle detection also depends on the size of the space-time image, and on the spatial ( x) and temporal ( t) resolution. Therefore, it is important to consider the limitations imposed by these factors on the precision of the angle detection. A space-time image of width w pixels and height h pixels (representing the space and time dimensions of the image respectively) would cover the distance x and time t as follows, The finest detectable change in the angle using the Radon transform for such an image is dependent on the number of streaks (n s ) in the image and the finest detectable change in the angle from a single streak (δ s1 ). Offsetting a single RBC streak trajectory by one pixel in the longer dimension of the streak can determine δ s1 (Figure 5). Assuming that only one streak in an image has changed the angle that is detectable with the Radon transform while other streaks maintained the same angle, the theoretically achievable finest detectable change in the angle (δ n ) would be, In other words, the finest detectable change in the angle improves with an increase in the number of streaks in a given image. For a given RBC streak spanning w s pixels in the spacedimension and h s pixels in the time-dimension, δ s1 and n s can be calculated as, where n s is the maximum number of streaks that can fit in the image. d s and d t are the minimum possible distance between two streaks in the space (in microns) and time (in seconds) dimensions, respectively. The streak angle (θ 0 ) can be calculated by application of the Radon transform on the image. With a known θ 0 , estimation of w s and h s for an image of width w pixels and height h pixels can be made using, The precision of angle detection is dependent on image dimensions (w and h), the detected RBC streak angle from the image (θ 0 ) and the number of RBC streaks that can fit in a recorded linescan (n s , which is in turn dependent on x, w, and d s ). The actual distance between RBC streaks can be greater than d s , making the actual number of streaks less than n s . Therefore, the realistic angle precision from a given image can be lower than δ n but higher than δ s1 .

Detection of change in the velocity is dependent on streak angle
Change in the blood velocity in a given vessel depends on a number of factors such as vessel size, vessel type (artery, capillary, vein), anatomical location, and stimulus type and can range from 0 to 80% of the baseline blood velocity (Drew et al., 2011;Shen et al., 2012). Therefore, the sensitivity of the algorithm needs to accommodate the expected change in the blood velocity without being computationally expensive. This can be achieved by optimizing the angle step-size used for the Radon transform in order to detect the desirable change in velocity. The fractional change in velocity ( v/v) can be defined as, Where v 1 is new velocity and v 0 is initial velocity. The corresponding change in RBC streak angle, δ = |θ 1 − θ 0 |, can be computed by combining Equations 9 and 15 and rearranging the terms, Thus, for the given RBC streak angle (θ 0 ), δ can be determined instantaneously because it is dependent on the desired v/v.
Likewise, the precision of v/v can be calculated from the value of δ used in the Radon transform and θ 0 . Note that the value of δ that can provide a given v/v is the highest for θ 0 = 45 • . Therefore, it is advantageous to perform line-scan imaging that yields ∼45 • streak angles. The smallest detectable δ, and therefore v/v, is dependent on the pixel dimensions of the image and θ 0 . δ smaller than the theoretically achievable finest detectable change in the angle (δ n ) cannot be resolved (see Equations 11-14 and Figure 6).

Application of the iterative Radon transform algorithm
The size of angle steps used in the calculation of Radon transforms determines the precision of angle measurements. This creates a trade-off between the achieved precision of angle measurements and the computational load. In the first implementation of the Radon transform to calculate streak angle (Drew et al., 2010), this trade-off was addressed by a two-step adaptive algorithm. First, relatively coarse angle steps of 1 • were used to achieve the approximate angle. Second, finer angle steps of 0.25 • around the measured angle were used. We refine the application of Radon transforms in a multi-step iterative manner to further optimize the balance between computational speed and precision of angle measurements.
The iterative Radon transform is started by first applying Radon transforms in sparse angle space and then finding the angle that gives the maximum variance in the image. In subsequent iterations, sets of angles are used that are less sparse than the preceding iteration and centered on the angle that yielded the maximum variance in the previous iterations (Figure 7). The smallest angle step-size (δ) to be used for a given image depends on the δ n , θ 0 , and the desired v/v (Equations 11,17). Since the variance at 0 • (vertical) is 0 when using a vertical Sobel filter, an angle selection of 0 • should be avoided.
The first iteration of the algorithm uses four Radon transforms with a 45 • angle step-size in order to span the full range of 180 • . In the second iteration, again four Radon transforms with a 45 • angle step-size are used with the center on the angle that yielded the maximum variance in the first iteration. Thus, the first two iterations give an angle resolution of 22.5 • . In the subsequent iterations, four Radon transforms with an angle step-size half that of the previous iteration are used. i.e., 22.5 • in the third iteration, 11.25 • in the fourth iteration, and so on. The number of iterations, i irt , required for achieving a given angle precision can be described as, Conversely, the precision achieved by a given number of iterations would be, Since each iteration performs four Radon transforms, the total number of Radon transforms (n irt ) used in the iterative Radon transform algorithm is, In order to achieve a given v/v sensitivity, δ values should be below the respective v/v curve, but above the δ n curve as determined by the image pixel dimensions. Note that the tolerance in the angle step-size values is highest at 45 • streak angle. For a given image size an angle step-size value below the solid line is not resolvable. δ n calculations are based on a space-time image with 1.19 μm/pixel spatial resolution and inter-streak distance of 4 μm (∼3.3 pixels). Open arrowheads represent the angle values with the maximum variance for given iterations, but smaller than the one from the previous iterations, and therefore were not carried forward to the next iteration. (C) Algorithm outputs for the first five iterations showing the angle range, angle step-size, and final angle output at the end of the iteration. The angles with the maximum variance for a given iteration are marked in bold. Angle value 47.5 • was found to have maximum variance in iteration 2, and was carried forward in proceeding iterations 3 and 4. This was because the maximum variances for the angle values in those iterations were smaller than that for

• (Compare variances represented by open and filled arrowheads in (B).
In this example, the iterative Radon transform offered a 12-fold increase in the angle precision and was still 4.5-times faster than 180 successive Radon transforms with angles spaced by 1 • (see Materials and Methods).

EFFECTS OF SOBEL FILTERING ON MEASURED STREAK ANGLE ACCURACY
Since for blood velocity the feature of interest in a given linescan image is the streaks resulting from plasma-RBC "junctions," it is important to ignore all the other elements of the image by minimizing them and/or enhancing the plasma-RBC junction. The vertical Sobel filter achieves this by removing low frequency changes in pixel luminance in the time dimension of the image.
Removing the aggregate mean of the image by subtracting it from each pixel of the image (whole-image demeaning) and removing the mean value across time at each spatial point of the line (temporal demeaning, vertical in our case) has been proposed previously for pre-processing of the image before performing Radon transforms (Drew et al., 2010). However, these methods do not effectively minimize the column-to-column variations in the luminance. Such differences are more pronounced when the image segments used for the Radon transform are small and/or the widths of the RBC streaks are wide for the given dimensions of the image segments ( Figure 8A). Although subtle, these ). (D) Sobel filtered images offer more accurate angle measurements when compared with vertically demeaned images. Mean ± SD of absolute difference in measured angles on simulated images are shown. Each angle value consists of results from n = 36 images with streaks at different spatial locations. Statistically significant differences are represented by a red asterisk for a given angle (p < 0.05, two-sample t-test).
column-to-column variations in luminance ( Figure 8B) lead to sub-optimal detection of streak angles. Application of Sobel filtering minimizes these variations and detected angles become more accurate ( Figure 8C) ). To further confirm the accuracy offered by Sobel filtering, we compared the differences in the measured angles from the actual angles on simulated line-scan images that were Sobel filtered versus vertically demeaned. A range of angles was tested with each angle having a sample size of n = 36 images (112 × 215 pixels, w × h), differing slightly in the spatial position of the streak. We found that angles measured on Sobel filtered images are more accurate when compared with the same image sequences that were vertically demeaned ( Figure 8D). For much longer timedimension images (200 × 1000 pixels, w × h), the improved accuracy on Sobel filtered images was retained, even though overall error rates were reduced with both types of filtering (data not shown).
In line-scan images acquired for blood velocity measurements, the Sobel filter acts as a high-pass filter by enhancing RBC streak margins while suppressing small line-to-line changes in pixel luminance. As a result, only the angle of RBC streak margins (plasma-RBC junctions) provides high variance in the Radon transform space. This improves the accuracy and precision of measuring the angle, and in turn velocity. Sobel filtering can be applied to line-scan images which have both thick (Figure 9) and thin (Figure 10) RBC streaks. In linescan images contaminated by brain movements that are due to transmitted pulsations from heartbeats (Figure 9B), vertical demeaning fails to suppress time-varying artifacts ( Figure 9C) and leads to incorrect RBC streak angle detection ( Figure 9E).
Apart from suppressing artifacts, Sobel filtering also enhances the RBC streak edges (Figure 9D) leading to an angle measurement with the iterative Radon transform algorithm that matches the actual RBC streak angle and therefore velocity ( Figure 9F). We also compared the angles measured after Sobel filtering vs. vertical demeaning over a long sequence of 35,000 contiguous line-scans (see Figure 9G). We found that ∼37% (n = 515) of image segments after vertical demeaning failed to accurately measure RBC streak angles, but no such errors were obtained when using Sobel filtered image segments ( Figure 9G).
Sobel filtering preserves the streak angle information in the images with thin, relatively horizontal RBC streaks (Figure 10). In line-scan images with very thin RBC streaks (Figure 10B), the Sobel filter enhances the RBC streak edges ( Figure 10D) when compared with vertical demeaning (Figure 10C), leading to detected angle values that are different from the vertically demeaned image (Figure 10E). The ∼0.6 • difference in the computed angle between the two filtering methods at ∼0.45 s in Figure 10E at the RBC streak angle of ∼83 • translates into a v/v of ∼8% (Equation 16), or ∼0.16 mm/s in this case ( Figure 10F). As the streaks become more horizontal, the difference between true and measured angle becomes similar after both vertical demeaning and Sobel filtering (see Figure 8D). Therefore, in the absence of time-varying artifacts, Sobel filtering provides little enhancement for images with very horizontal streaks.

IMPROVED SPEED AND PRECISION OF VELOCITY MEASUREMENTS WITH ITERATIVE RADON TRANSFORMS
Blood velocity measurements can only be as precise as the angle step-size used in the Radon transform. But decreasing the stepsize to improve precision results in additional Radon transforms and increased computation time. The iterative Radon transform method helps overcome this trade-off. The total number of Radon transforms required by the traditional non-iterative method (n trt ) spanning the full range of 180 • for the step-size (δ) would be, Compared to the traditional non-iterative method, the method described here requires an order of magnitude fewer Radon transforms. Based on Equations 18, 20, and 21, the relationship between the total number of Radon transforms required by the traditional non-iterative method (n trt ) and the iterative method (n irt ) is, The angle precision (δ) of 1 • requires 180 Radon transforms with the non-iterative method but only 28 using the iterative method. Similarly, an angle precision of 0.01 • would require 18,000 Radon transforms if performed non-iteratively, but can be resolved with only 56 Radon transforms using the iterative method (Figure 11). The iterative Radon transform method overcomes the speedprecision trade-off by starting with sparse sampling of the angles and then sampling progressively finer angles in the subsequent iterations centered around the angle with the highest variance (see Materials and Methods).

DISCUSSION
The hybrid algorithm presented here performs pre-processing of the image followed by iterative application of Radon transforms. Pre-processing of the images with Sobel filtering has an advantage over vertical demeaning in that it removes the slow time-varying artifacts in addition to time-invariant artifacts. This is especially useful in large animal preparations where brain movement artifacts due to respiration and heartbeats are sometimes difficult to eliminate. Such artifacts lead to line-scan images with discontinuous RBC streaks and/or uneven or varying background luminance. Additionally, Sobel filtering also enhances the RBC streak margins. Other types of filtering using operators similar to Sobel, e.g., Prewitt or Scharr operators can also be used with equivalent results. Sobel filtering results in an increased contrast at the plasma-RBC junctions. Because the Sobel operator we use is oriented vertically, the margins of relatively horizontal streaks are enhanced more than vertical streaks. Completely vertical RBC streaks appear in cases of stagnant blood flow (Kleinfeld et al., 1998). In these situations, any vertical operator (including Sobel, vertical demeaning) will suppress the streaks and therefore should be avoided. The detection of the correct angle using the iterative Radon transform is dependent on the SNR of the Sobel filtered image. We found in our simulations that almost vertical stripes (angle < 1 • ) can be efficiently detected using the iterative Radon transform algorithm after Sobel filtering. However, the detection FIGURE 11 | Iterative Radon transform method markedly reduces the required number of Radon transforms (n irt ) to achieve a given angle precision (δ) when compared with non-iterative method (n trt ). The number of Radon transforms in the iterative algorithm increases logarithmically as δ becomes finer, compared to a linear increase for the traditional, non-iterative counterpart. As a result, there is an exponential increase in the computation speed for the iterative algorithm compared to the traditional non-iterative Radon transform methods. Results are based on Equations 18-22.
is less robust to noise than in the case of stripes with higher angle values. More generally, velocity measurements using the algorithm described here will produce the best results provided that the angles of the RBC streaks are close to 45 • (see Materials and Methods).
Changes in velocity due to heart rate pulsations can be substantial (see Figures 9, 10) and may exceed the changes due to sensory stimulus-evoked responses. However, the v/v induced by heart rate pulsations has the frequency range of 3-8 Hz which is much higher than a typical <1 Hz v/v evoked by a given sensory stimulus (Kleinfeld et al., 1998;Chaigneau et al., 2003;Berwick et al., 2008;Drew et al., 2011;Shen et al., 2012). Therefore, such physiological v/v can be band-pass filtered using Fourier transform methods to generate more accurate measurements of stimulus-evoked v/v. A low-pass filtering can also be achieved by averaging the velocities over time-windows that span multiple cardiac cycles. Alternatively, the collected line-scan data can be consistently phase-locked to the cardiac cycle by triggering line-scan imaging at a particular phase of the cardiac cycle, e.g., using QRS complex as a trigger. The angle step-size (δ) we used for Figures 9, 10 were ∼0.01 • and ∼0.005 • , respectively, in order to achieve the desired detection sensitivities of 0.1% v/v (see Equation 17). However, if the experimental paradigm does not demand detection of minute changes and/or wider range in the velocities, the user may choose a bigger δ and/or narrow starting angle range for even faster computation speeds.
Iterative application of the Radon transform resulted in an order of magnitude faster determination of blood velocity when compared with other algorithms of blood velocity measurements.
When testing simulated line-scan data with RBC streak shifts ranging from 10 to 80 pixels per line and using ∼0.1 • angle precisions on a 200 × 15,000 image, our algorithm was ∼6.5× faster than the correlation-based method (Kim et al., 2012). However, the correlation-based method offers superior velocity measurements in some instances of near-horizontal RBC streaks. A spatial frequency-based analysis claims 1.25× faster performance over Radon transform when using 1 • angle stepsize spanning 180 • (Autio et al., 2011). Our algorithm achieves this precision in only 28 Radon transforms (Equation 22), making it 6.4× faster than the Radon transform when applied in a non-iterative manner. Thus, the advantage offered by the iterative implementation makes our algorithm 5.14× faster compared with a spatial frequency-based analysis. The "fasttrack" algorithm that uses global energy minimization technique claims superior velocity measurements than Radon transform in near-horizontal RBC streaks, but is generally 5× slower than traditional Radon transform (Deneux et al., 2011(Deneux et al., , 2012. Conservatively estimating a 6.4× speed improvement offered by the iterative implementation of Radon transforms over its non-iterative counterpart makes our algorithm 32× faster than the "fasttrack" algorithm. Furthermore, Sobel filtering followed by very fine angle measurements in an iterative manner eliminates the likelihood of detecting false angles that can otherwise occur when line-to-line discontinuities of RBCs are present in space-time images. The global energy minimization technique appears to have a comparable accuracy to the standard Radon transform method when RBC streaks are not near horizontal (Deneux et al., 2011). The apparently poor performance of the Radon transform algorithm on detecting near-horizontal RBC streaks (high velocities) in the study by Deneux and colleagues may be due the use of relatively coarse angle steps and incomplete removal of slow time-varying luminance variations in the images that are unrelated to the flow of RBCs.
In summary, we present a method of accurately and precisely measuring blood velocities from line-scan images in a computationally efficient manner. First, Sobel filtering is used to remove slow time-varying artifacts and enhance the edges of RBC streaks. This leads to accurate estimation of RBC streak angles and thus blood velocity measurements. Second, the iterative Radon transform minimizes the required number of Radon transforms needed to precisely estimate the RBC streak angle. This leads to fast RBC streak angle measurements. The near-ideal precision attained with the iterative Radon method together with the high accuracy provided by Sobel filtering ensures that relatively small differences in baseline vs. stimulus-evoked RBC streak angles can be measured. This is a non-trivial advance because data collected from in vivo experiments include movement artifacts. Such data may otherwise be discarded, or if used would lead to erroneous measurements. More generally, the individual components of our algorithm can be used independently for a broad range of particle velocity applications. For example, Sobel filtering of noisy spacetime images collected for tracking the movement of any particle can lead to improvement in the SNR by enhancing the particle streak edges. The concept of iterative image analysis can be used to speed-up other algorithms that have been used to calculate the slope of RBC streaks (Drew et al., 2010;Duncan et al., 2010, Autio et al., 2011. Our hybrid algorithm (in its entirety or as individual components) can also be used for real-time implementation of velocity calculations in high-speed imaging applications that use resonant or acousto-optic deflector scanners.