1. Introduction

Optical coherence tomography (OCT) is a 3D biomedical imaging technique based on optical interference to achieve depth-revolved mapping with a penetration depth on the order of 1-2 millimeters [1–3]. Through a process called coherence gating, OCT analyzes the interference between light reflected from a specimen and a reference beam to determine the optical properties of a sample in a depth-resolved manner.

Doppler OCT is a special application of this technique. It uses the time-varying phase of the reflected light from the specimen to determine the velocity component at every voxel in the direction of the probing OCT beam. The time-varying signal obtained from Doppler OCT can be represented on the complex plane. An idealized noise-free signal rotates about the origin and thus its phase angle linearly increases in time. Doppler OCT uses the rate of this signal rotation to determine velocity, as the rotation rate is proportional to the velocity. Finding the rate of rotation in the ideal signal is simple; one can determine the change in phase or angle that each data point makes about the origin per unit time, often in conjugation with an autocorrelation to reduce noise [4–6]; one can determine the mean or maximum frequency of the signal after a frequency decomposition [7,8].

In real signals, however, the center of rotation (COR) of the signal is not always the origin. This will occur in a specimen when the voxel under consideration contains both static and mobile scattering components. The static scattering component introduces a superimposed offset in the signal [9–11] that leads to velocity underestimation and reduced estimator performance [12,13]. Some examples from Doppler OCT data of rodent cerebral blood flow clearly show this complication in real signals (Fig. 1).

 

Fig. 1 Examples of non-origin CORs in real Doppler OCT signals. Multiple rotations in the top row correspond to higher velocities, while partial rotations in the bottom row represent lower velocities. All signals are shown on the same scale (a.u.). Re, real axis; Im, imaginary axis.

Download Full Size | PDF

High-pass filtering has been widely used to determine the COR as the mean of the signal data points [10,14–16]. This method works well in cases of large signal rotations (i.e., high velocities) but is biased for smaller rotations (Fig. 2). This bias makes the COR-corrected signal seem to span a wider angle and thus have a faster rate of phase change than the true values, leading to potential overestimation of velocity [17]. Maximum and mean frequency-based methods are also susceptible to static clutter [8], and the latter case at least demonstrates bias when high pass filtering is applied [8,18]. Other methods for clutter rejection also show some limitations [15,16,19–21]; examples include a second-order autoregressive model which has non-ideal performance for large clutter signals and can demonstrate spectral leakage of the clutter component [15,22], and the mean pursuit method which, despite its high computational cost, may still remove non-clutter components or fail to remove some clutter components [16]. The few data points available in the signal pose what is often regarded as one of the most difficult challenges for designing clutter filters [23].

 

Fig. 2 Illustration of the bias in the high-pass filtering method to determine CORs. Re, real axis; Im, imaginary axis.

Download Full Size | PDF

Here, we propose two novel methods to more accurately determine the COR and thereby increase the accuracy and dynamic range of Doppler OCT measurement, even when the velocity is so low that the signal makes only partial rotations. First, the COR can be determined as a point on the complex plane such that the variance in the distances from the determined COR to the data points is minimized. This method, termed variance minimization, will correctly give CORs for both large and small signal rotations.

However, even variance minimization may be biased when noise is large enough to move the variance minimizing point toward the data points and away from the true COR. Therefore, we propose and test a second method, termed offset extrapolation, based on the statistical characterization of this potential bias and its application to extrapolate the true COR from the high-pass filtering and variance minimization-determined CORs.

This paper describes the implementation of these two algorithms, numerical simulations to test and compare their performance with the current standard method of high-pass filtering, specifically comparing the accuracy of both COR determination and velocity measurement, and applications to experimental data to demonstrate the effects of these novel algorithms.

2. Methods

2.1 Numerical data

We used 16 time points (i.e., 16 consecutive A-scans) for our numerical simulations. We chose a set rotation for the signal (larger rotations corresponding to larger velocities). For characterization of the bias, we tested 40 different rotations from 0 to 2π, and fitted a bias or offset function to the data. For validation, we specifically tested π, 5π/4, 3π/2, 7π/4 and 2π rotations, which corresponds to flow velocities of 5.1-10.2 mm/s for our OCT system with 1,310-nm center wavelength and 147,000 A-scan/s speed. We then added Gaussian white noise to each data point independently. In implementing the offset extrapolation method, we used 8 different noise levels, with standard deviations of the noise being 1/8, 2/8, 3/8, 4/8, 5/8, 6/8, 7/8, and 8/8 of the radius of rotation (ROR) of the OCT signal. The ROR was set to 1 in our numerical simulation. For the performance evaluation, we tested three noise levels: 0, 1/16, and 1/8 of the ROR. Our tests were independent of the orientation and offset in the simulated signal, since the COR algorithms always determined the COR point with respect to the signal, independent of orientation and offset. For implementing the offset extrapolation method, we tested 1000 data point for each rotation and noise case. We tested the COR algorithms with each 500 generated signals for the 1/16 ROR noise case, and 1000 signals for the 1/8 ROR noise case for each of the different rotations.

2.2 Spectral-domain OCT system

A commercial spectral-domain OCT system (Thorlabs, Newton, NJ, USA) was used for experimental data acquisition. The system uses a large-bandwidth near-infrared (NIR) light source with a center wavelength of 1310 nm and a wavelength bandwidth of 170 nm, which leads to high axial resolution (3.5 μm). The 5x NIR objective lens (Mitutoyo) was used for lateral resolution of 7 μm. The focal plane was located ~120 μm below the cortical surface. The system uses a high-speed 2048-pixel line-scan camera to achieve 147,000 A-scan/s with a relatively large imaging depth. Thus, the temporal sampling rate of Doppler OCT signal was 6.8 μs. The field of view (FOV) had 512 (X) x 512 (Y) x 1024 (Z) voxels with the transverse and axial sampling rates of 3 μm and 3.5 μm, respectively, leading to 1.54 mm x 1.54 mm x 3.58 mm imaging volume. We repeated 8 A-scans at each (x,y) position, leading to the temporal signal duration of 54 μs.

2.3 Animal preparation and imaging

Five mice were used for this study. On each mouse, craniotomy surgery was performed at 9 weeks of age. The animal was sedated with 3% isoflurane and sedation maintained with 1.5% isoflurane, delivered with 1 L/min O₂. An approximately 10 mm sagittal incision was made on the surface of the head, terminating at least 3mm behind the eyes. The incision was symmetrically laterally expanded to expose a 10mm-diameter surface of the dorsal skull. A custom metal head post was cemented to the cortex. Once cured, a circular hole of diameter 3 mm was created in the right parietal bone. A glass aperture consisting of one 5 mm x 0.1 mm and two 3 x 0.1 mm circular cover glasses were affixed, leaving a sealed 3mm aperture which allowed direct longitudinal optical viewing of the cortex.

After surgery, the animal was placed into the OCT system for Doppler OCT imaging of cerebral blood flow. The animal was kept under isoflurane anesthesia during imaging. Oxygen saturation, pulse rate, and temperature were continuously monitored with pulse oximetry and a rectal probe throughout the surgical procedure and experiment. The body temperature was maintained at 37°C and the pulse rate remained within the normal range of 250–350 pulse/min.

All animal-based experimental procedures were reviewed and approved by the Institutional Animal Care and Use Committee of Brown University, according to the guidelines and policies of office of laboratory animal welfare and public health service, National Institutes of Health.

3. Results

3.1 Variance minimization method

In theory, on the complex plane, the COR should have equal distances from Doppler OCT signal data points. Thus, in practice, we can determine the COR as a point where the variance in the distances is minimized [11]. We derived an equation to find the point in the complex plane ($A+Bi$) for a given set of Doppler OCT signal data points, $R\left(t_k\right)=a_k+b_{k}i$.

 

Fig. 3 Examples of CORs determined by the high-pass filtering and variance minimization methods. Re, real axis; Im, imaginary axis.

Download Full Size | PDF

In contrast, when noise is very large, there is a chance for the variance minimization method to be biased as illustrated in Fig. 3c. The illustrated example shows that the COR determined by Eq. (1) is biased toward the data points from the true COR. As noise increases, the variance minimization-determined point moves closer to the mean of the data points (i.e., the high-pass filtering-determined COR). In summary, the CORs determined by the high-pass filtering and variance minimization methods follow a pattern: Firstly, reducing the rotation increasingly biases the high-pass filtering COR away from the true COR and towards the signal. Noise does not affect this bias. Secondly, small rotations and large noises may bias the variance minimization COR away from the true COR towards the signal, but less than that of the high-pass filtering COR. The patterns inherent in these potential biases have motivated us to develop another method as follows.

3.2 Offset extrapolation method: determining the offset function

First, we systematically characterized how the above biases in the high-pass filtering and variance minimization methods depend on the degree of rotation and noise level, using the numerical data with known true values of the COR, rotation or velocity, and noise level. This characterized relationship will serve as a lookup table in our second method, termed offset extrapolation. While the lookup table receives the rotation and noise level as its input, the input values do not need to be precise because the offset extrapolation method will use an iterative search to find their accurate values as described in Section 3.3.

The noise level was defined by the ratio of the standard deviation of added Gaussian noise to the radius of rotation (i.e., overall amplitude of the OCT signal’s mobile component). For each rotation and SNR, Doppler OCT signals have been generated as described in Section 2.1 (n = 1000 per rotation and noise case with individual random noise generation). From each signal, we determined two CORs by using the high-pass filtering and variance minimization methods; then, the CORs were ensemble-averaged (over n = 1000 per rotation and noise case). Finally, we quantified the offset between the two CORs and the true COR as shown in Fig. 4(a). The offset was obtained for each rotation and SNR, resulting in the offset as a function of two parameters (Fig. 4(b)).

 

Fig. 4 (a) Definition of the offset between the high-pass filtering-determined COR (green), variance minimization-determined COR (blue), and the true COR (black). Re, real axis; Im, imaginary axis. (b) The characterized offset as a function of the degree of rotation and noise level. The rotation represents flow velocity when the OCT scan speed and data point number are fixed. The noise level is presented in the standard deviation (SD) of added random noise divided by the radius of rotation (ROR).

Download Full Size | PDF

The offset function takes the rotation of the signal and the noise level as inputs and returns the offset. This function is shown to have an approximately quadratic relationship to noise (increasing with larger noise) and a reciprocal relationship with rotation (increasing with smaller rotations). Note that the obtained offset function or the lookup table only depends on the geometry of a given Doppler signal, not on any system parameter (like A-scan/s); thus, it is universally applicable for different systems.

3.3 Offset extrapolation method: iterative search

When the true COR is known, the degree of rotation in a given Doppler OCT signal can be determined by finding the total angle spanned by the signal from the COR. The radius of rotation (ROR) can be estimated as the mean distance to the data points from the COR and a variance in the distance computed. Dividing the variance by the mean distance will give an estimate of the noise level as defined above. These rotation and noise levels measured with respect to the true COR will serve as an accurate input for the offset function which in turn will enable us to find the true COR again. However, in practice, all of the true COR, true rotation, and true noise level are unknown. Therefore, we have implemented an iterative search so that we can find the ‘true’ COR from a provisional COR, without any priori information.

As illustrated in Fig. 5, for a given signal, we determine two CORs by the high-pass filtering and variance minimization methods, choose a provisional COR, estimate the rotation and noise level with respect to the provisional COR, input them to the offset function, and use its result and the two CORs to extrapolate a ‘true’ COR. Then, this process is iterated until the discrepancy between the old and new ‘true’ CORs becomes smaller than a pre-defined threshold. The threshold we used was determined independently for each voxel and was taken to be 1/1000 of the maximum distance between the data points and their mean.

 

Fig. 5 Diagram of the iterative process used to determine the offset extrapolation method. HP, high-pass filtering; VM, variance minimization.

Download Full Size | PDF

Our implemented iterative algorithm has demonstrated convergence in the range of parameters tested (see Section 2.1. for detailed parameter range). When a provisional COR is chosen to be very far from the origin and the data points, it makes the discrepancy between the old and new CORs increase indefinitely, as the spanned angle and the noise level follow an inverse relationship with an increase in the discrepancy. However, even in such cases, the iterative search does not diverge as our offset function increases faster with respect to a decrease in the rotation than it decreases with respect to a decrease in the noise level.

An example of Fig. 6 shows that this method determines the COR more accurately than both the high-pass filtering and variance minimization methods, especially when the signal data exhibits a partial rotation and large noise. A more systematic evaluation of the method is presented in later sections.

 

Fig. 6 Examples of four CORs for different rotations and noise levels: true, high-pass filtering (HP)-determined, variance minimization-determined, and offset extrapolation-determined CORs. Re, real axis; Im, imaginary axis.

Download Full Size | PDF

Both methods were implemented in MATLAB, and our software codes have been uploaded to Brown Digital Repository with the DOI of 10.7301/Z0K35S6G for public access.

3.4 Comparative algorithm performance

We compared three methods (high-pass filtering, variance minimization, and offset extrapolation) in terms of their accuracy in determining CORs. The accuracy was measured as an absolute percentage of the ROR (the overall amplitude of the OCT signal’s mobile component). The accuracy was tested for various degrees of rotation (or flow velocities) and noise levels. As presented in Fig. 7, both the variance minimization and offset extrapolation methods result in higher accuracy than the previously used high-pass filtering method, by up to 60 percent points. As expected, this improvement is more significant with partial rotations (low velocities). Using the OCT scanning speed of 147,000 A-scan/s and 8 A-scans per (x,y) position, the partial rotations of $\pi$ -$2\pi$ correspond to axial flow velocities of 5.1 – 10.2 mm/s, which are within the known physiological blood flow speed range [5,24,25]. When noise is as large as 10% of the ROR, the accuracy of the variance minimization method drops to ~90% in the partial rotations, but the offset extrapolation method maintains near 100% accuracy.

 

Fig. 7 Accuracy in COR determination. The example velocity values correspond to the rotation values when 147,000 Ascan/s and 8 Ascan/position are used. SD, standard deviation of simulated random noise; ROR, radius of rotation.

Download Full Size | PDF

The performance of the three methods was also tested in terms of the accuracy of velocity determination, in conjugation with different Doppler analysis algorithms: the Kasai [4,25], the autocorrelation function fitting, and the maximum likelihood estimation (MLE). In short, the Kasai algorithm focuses on the phase change in the autocorrelation function at the first time lag, the autocorrelation function fitting algorithm performs linear fitting on multiple phase data of the autocorrelation function to determine the rotation rate (specifically up a lag of half of the total time), and the MLE finds the maximum amplitude frequency in the signal power spectrum with no autocorrelation. As can be seen in Fig. 8, compared to the high-pass filtering method, both the variance minimization and offset extrapolation methods result in up to 60% higher accuracy, although the variance minimization method shows ~5% lower accuracy than the offset extrapolation in large-noise signals. The Kasai method is least affected by the COR determination accuracy, so we used it in the following in vivo data analysis.

 

Fig. 8 Accuracy in velocity measurement. The example velocity values correspond to the rotation values when 147,000 Ascan/s and 8 Ascan/position are used. SD, the standard deviation of simulated random noise; ROR, the radius of rotation.

Download Full Size | PDF

3.5 Application to in vivo Doppler OCT imaging

We tested how the proposed clutter rejection methods work on in vivo data acquired from the rodent cerebral cortex. First, compared to the high-pass filtering method, the proposed methods are expected to be more robust against noise as shown in Fig. 8. To visually confirm this, we compared en face Doppler OCT velocity maps at various depths between the high-pass filtering and variance minimization methods. As can be seen in the example of Fig. 9 (dashed white boxes), the proposed method results in much less noise, especially in the deeper cortex (Fig. 9(c)) where the Doppler OCT signals are believed to have larger noise due to the multiple scattering and overall signal attenuation. In addition, the proposed methods are expected to produce more accurate results especially for low velocities (i.e., smaller rotations as shown in Fig. 8), which in turn leads to better visualization of small blood vessels in Doppler OCT velocity maps. The example of Fig. 9 also qualitatively shows this tendency; small blood vessels like arterioles (white arrow in (a)) or even some capillaries (white arrow in (b)) were visible with the variance minimization method.

 

Fig. 9 Examples of en face OCT angiogram and Doppler OCT velocity maps. Every en face map was obtained by either maximum intensity projection (angiogram) or mean intensity projection (Doppler) over ± 35 μm depth. Scale bar, 100 μm.

Download Full Size | PDF

To evaluate the clutter rejection methods in a more quantitative manner, we tested how the absolute blood flow is conserved within the same blood vessels, because it is difficult to obtain the exact values for blood flow velocities in such animal data. The absolute blood flow can be obtained in the unit of μL/min via an area integral from an en face slice of the 3D Doppler OCT velocity map [26]. The flow should be conserved within the same vessel segment. Flow conservation checks were carried out by comparing blood flow values at two different locations (depths) of a penetrating artery or draining vein, for 5 veins and 5 arteries in each of the 5 mice (50 vessel segments in total). Different algorithms result in different Doppler OCT velocity maps and thus different flow values, as shown in the example of Fig. 10.

 

Fig. 10 Examples for en face Doppler OCT velocity maps and blood flow values, for different clutter rejection methods. Scale bar, 100 μm.

Download Full Size | PDF

From the two flow values measured at different locations, the conservation error was quantified by {(larger value) – (smaller value)}/(larger value). Statistically, no clutter rejection produces the highest error for flow conservation (Fig. 11(a)). The statistical result also suggests that the offset extrapolation and variance minimization methods outperform the currently-used high-pass filtering method. Although the offset extrapolation method was slightly more accurate than the variance minimization in the numerical simulation (Fig. 8), they produce statistically insignificant differences in performance when applied to the animal data.

 

Fig. 11 Statistics of the flow conservation error (a) and the absolute flow (b). *p<0.05, ***p<0.001 (paired t-test). Data are presented as mean ± SE. NS, Not Significant.

Download Full Size | PDF

Figure 11(b) shows the underestimation of absolute blood that occurs when no clutter rejection filter is used. We can also see the overestimation present in the high-pass filtering and, to a lesser extent, variance minimization. Overestimation is less pronounced for the offset extrapolation method. However, as the exact true values are unknown in the in vivo data, it is difficult to quantify the degree of these biases or to know whether the offset extrapolation method produces an overestimation or underestimation of the flow velocities.

4. Discussion and conclusion

This study describes two novel methods for Doppler OCT clutter rejection and demonstrates the improved accuracy in both numerical tests and in vivo imaging. As can be seen in the numerical test (Figs. 7 and 8), the proposed methods are specifically accurate in lower velocities (partial rotations; corresponding to lower than 7.6 mm/s in our OCT system), improving both COR and velocity accuracies by up to 60%. Higher accuracy may become increasingly important with recent advances in OCT speed, as faster OCT scanning will lead to smaller rotations in the Doppler signal when the other parameters are identical. The improved accuracy in low velocities (i.e., extended dynamic range) may allow for more robust detection of peripheral blood vessels whose flow velocities tend to be lower, and thus help avoid the blood vessel size underestimation that is currently observed [12]. In addition, when noise is as large as 10% of the ROR (which is common in real in vivo Doppler OCT data), the offset extrapolation always, though slightly, outperforms the variance minimization method in the numerical tests but not significantly with the in vivo imaging data (Fig. 11(a)), perhaps due to the higher variance of the offset extrapolation estimator leading to a greater overall MSE. Finally, while no clutter rejection leads to velocity underestimation and the high-pass filtering method tends to lead to overestimation, our methods produced intermediate flow values (Fig. 11(b)), which might be associated with the improved accuracy.

With a longer acquisition time, the complex OCT signal will tend to rotate more than 2π (a full cycle) so that one can expect the high-pass filtering to work better. However, “longer” here can be a relative concept. More accurately speaking, a longer acquisition time will lower the bottom limit of the velocity range for which the high-pass filtering works accurately, but the HP filtering still would not work well for lower velocities than the new limit. Meanwhile, a shorter acquisition time is generally advantageous in many aspects; for example, it allows for faster dynamic imaging (i.e., more volume data) with repeating scans to trace temporal changes in CBF against time for functional studies. Thus, there is a tradeoff between acquisition time and dynamic range. It is noteworthy to mention that even with a longer acquisition time, there always exists a velocity range for which the variance minimization and offset extrapolation methods work more accurately than the high-pass filtering method.

Although the offset extrapolation method produces the most accurate results in all of our numerical tests, it is based on the iterative search which requires the calculation of the high-pass filtering and variance minimization CORs and thereby is the most computationally intensive of the three algorithms. Given this, along with the statistically insignificant difference in the in vivo imaging result, the variance minimization method may be a fair compromise between accuracy and computational intensity. If one wants to reduce the computational overhead of offset extrapolation, various ideas could be tested, such as fitting the mismatch of the provisional and extrapolated CORs to a quadratic function, and using the zeros of that function to produce a better estimate of the locations of the extrapolated CORs, instead of using the simple bisection method adopted here.

Offset extrapolation also showed the greatest variance in the determined COR in our numerical test. It might be attributed to biases that we observed in the numerically determined offset function, especially at low noises. A more accurate determination of the offset function may improve the precision for very small rotations. Also, the blood flow conservation test might not be ideal for in vivo determination of the accuracy. There is a chance for a method to falsely produce high flow conservation accuracy when the method has a systemic error (e.g., consistently overestimate the flow). However, since a recent study validated the overall accuracy of Doppler OCT measurements by utilizing the traditional hydrogen clearance method to measure cerebral blood flow [17], we do not expect that the proposed two methods for COR determination would specifically produce a systematic error.

Finally, of particular note was that the Kasai method, even when combined with the previous high-pass filtering method, resulted in relatively robust accuracy in our numerical test (Fig. 8), despite large inaccuracies in the COR (Fig. 7). In our numerical test, we first removed the COR from the signal and then obtained the autocorrelation function from the COR-corrected signal. It seems that, for very large biases in the determined COR, the autocorrelation function with a single time lag (used for Kasai) eliminates the large jump that occurs in the phase angle. This might help to mitigate the velocity overestimation otherwise present in the method. The variance minimization and offset extrapolation methods, however, are much more robust to the autocorrelation time lag chosen; for instance, they perform similarly at the unit lag as well as the half lag.

Even in conjunction with Kasai, the three methods exhibited a difference in their performance. Although an exact theoretical basis underlying the difference should be further studied, the large phase angle jump due to incorrect COR estimation in the high-pass filtering method might be mitigated when using an autocorrelation with unit lag. Nonetheless, the improvement by that mitigation was less pronounced, apparently due to the lower bias of the variance minimization and offset extrapolation methods in determining the true COR.

In conclusion, we have demonstrated new clutter rejection methods for Doppler OCT. While the old algorithms exhibit biases, our methods are designed to eliminate the source of this bias. This should make Doppler OCT robust to different circumstances where the static component makes a stronger or weaker contribution to the signal, an advantage over simple high-pass filtering of OCT signals. We expect our algorithms to enable Doppler OCT to present a truer picture of the absolute blood flow measurements made across different in vivo scenarios.