Real-time bulk-motion-correction free Doppler variance optical coherence tomography for choroidal capillary vasculature imaging.

In this paper, we analyze the retinal and choroidal blood vasculature in the posterior segment of the human eye with optimized color Doppler and Doppler variance optical coherence tomography. Depth-resolved structure, color Doppler and Doppler variance images are compared. Blood vessels down to the capillary level were detected and visualized with the optimized optical coherence color Doppler and Doppler variance method. For in-vivo imaging of human eyes, bulk-motion induced bulk phase must be identified and removed before using the color Doppler method. It was found that the Doppler variance method is not sensitive to bulk-motion and the method can be used without correcting the bulk-motion when the sample-movement-induced velocity changes gradually. Real-time processing and displaying of the structure and blood vessel images are very interesting and is demonstrated using a dual quad-core Central Processing Unit (CPU) workstation. High resolution images of choroidal capillary of the vasculature network with phased-resolved color Doppler and Doppler variance are shown.

1. Introduction

Optical coherence tomography (OCT) is a powerful interferometric technology used to obtain cross sectional images of tissue noninvasively with micrometer resolution, millimeter penetration depth and video-rate imaging speed [1]. OCT has become a valuable tool in a number of medical fields, especially in ophthalmology due to its non-contact, high resolution nature. Extension of the OCT technique to functionally image the human eye is of great interest. Doppler optical coherence tomography (DOCT) is one kind of functional extension of OCT which combines the Doppler principle with OCT and provides in-vivo imaging of blood vessels, blood flow direction, and flow speed [2-17].

Blood flow in humans is complex, and the flow is pulsatile. The flow speed in a blood vessel is not a constant value, and changes from the center of the blood vessel to the edge. A parabolic equation may describe the distribution of the steady blood flow along the blood vessel radius. This kind of steady flow is usually called laminar flow. Laminar flow may breakdown into a turbulent flow when the velocity becomes sufficiently high [18]. In practice, the Doppler frequency produced by a single target will produce a spectrum of Doppler frequencies or a series of frequency shifts instead of a single frequency. This spectrum broadening is attributed to several sources, such as the cone-geometrical focusing beam, Brownian motion and speckle. Brownian motion dominates the broadening of the Doppler spectrum at low flow speed, and probe-beam geometry dominates at high flow speed [19]. Depending on the Doppler frequency shift information obtained, we can display this information in a color Doppler method or a variance method. In the color Doppler imaging method, the average Doppler frequency shift is displayed as color images, and the negative and positive averaged frequency shifts are displayed in different colors [2-4], which provide quantitative information on the flow speed and flow direction. In a Doppler variance image, the variance or the standard deviation of the Doppler frequency shift is displayed [5], which can be used to quantify Brownian motion or measure transverse flow [19].

In early time domain OCT systems, the Doppler frequency shift was obtained via a spectrogram method which used the short time fast Fourier transformation (STFFT) or wavelet transformation [2-4]. However, the spectrogram method allowed detection of high speed flow only, and a phase-resolved method was introduced to solve this problem [6]. Phase-resolved color Doppler has been successfully used to analyze retinal blood flow in real-time [7,8]. Recently, several groups has successfully extended phase- resolved DOCT method to image blood vessel networks in the choroid [9-16]. Yasuno et al. demonstrated an intensity based binarization method- scattering optical coherence angiography (S-OCA) to indentify the blood vessels [9]. The authors demonstrated fine vasculature in the retinal and choroid layers with this S-OCA method [10]. Wang et al. proposed a method called optical microangiography (OMAG) technology to separate the static and moving signals with a modified Hilbert transform [10-14]. By stitching multiple small scanning areas, high-resolution choroid blood vessel images down to the capillary level were demonstrated with OMAG technology [12]. By applying the OMAG algorithm along the slow scanning axis, the authors achieved ultrahigh sensitivity imaging to image the slow flows within capillaries [13]. However, the algorithm can be applied only after removing bulk-motion using the phase-resolved method and the introduction of the additional Hilbert transform increases the algorithm’s complexity and computation time. Imaging multiple small areas will increase the imaging time, and registering several small area images into one image is technically challenging. Another method named Joint Spectral and Time Domain OCT also demonstrated the detection of capillary level blood vessels [15]. However, this method needed a high scanning density which may greatly increase imaging time. Fingler et al. proposed a phase variance method that uses a sliding window of 10 sequential B-scans to calculate intensity averages and phase variance for retinal blood vessel imaging [16]. Recently, our group has demonstrated imaging of retinal and choroid blood using phase-resolved Doppler variance imaging [17].

In this paper, we show that imaging the choroid capillary vasculature network is possible using phased-resolved color Doppler and Doppler variance methods without sacrificing imaging time and data processing time. For in-vivo applications, especially ophthalmology, axial sample movement will induce bulk-motion and change the detected Doppler frequency. The axial sample movement induced bulk-motion must be corrected before applying the optical Doppler tomography algorithm. Histogram-based methods are often used to determine the bulk-motion induced bulk phase. However, this method may also introduce additional artifacts, and it is computationally intensive. We show that Doppler variance is not sensitive to bulk-motion and the method can be used without correcting the bulk-motion when the sample-movement-induced velocity changes gradually. We also demonstrate real-time processing and displaying of structure and blood vessels with a dual quad-core workstation. High resolution imaging of the choroid capillary vasculature network with phased-resolved color Doppler and Doppler variance is shown.

2. Method

Phase-resolved Doppler OCT has become the prominent method to extract the average Doppler frequency, and is based on the detection of the phase difference between adjacent A-lines [6,18].

where is average Doppler frequency shift and is the power spectrum of the Doppler frequency shift, and are the phase for the signal at depth z of and A-line. Based on autocorrelation technology, Eq. (1) can be rewritten as [20]:

where arg(0) is the argument function, and are the real and imaginary parts of the complex data . This so-called autocorrelation algorithm has shown great success in the fields of both ultrasound Doppler imaging and Doppler OCT [6,20].

Doppler variance, an extension of the OCT/ODT technique uses the variance of the Doppler frequency spectrum to map the flow. Doppler variance has benefits of being less sensitive to the pulsatile nature of the blood flow, less sensitive to the incident angle, and may be used to obtain the transverse flow velocity [5,17,19]. If σ denotes the standard deviation of the Doppler spectrum, the Doppler variance can be obtained [5,20]:

With the help of autocorrelation technology, the variance can be expressed as [20]:

Usually the equations are used together with averaging to improve signal to noise ratio (SNR). Averaging can be performed in the lateral direction (temporal direction) so that Eq. (2) and Eq. (4) become [5,6,20]:

where J is the number of A-lines that are averaged. Averaging could also be performed in both the lateral and depth direction, and Eq. (2) and Eq. (4) become [21,22],

where J is the number of A-lines that are averaged and N is the number of depth points that are averaged and the values of J and N are application dependent. Generally, a larger J and N will increase SNR, but also increase computation time and decrease resolution.

3. Axial sample movement and its impact on the color Doppler and Doppler variance

For in vivo applications of Doppler OCT, axial sample movement will induce bulk-motion and change the detected Doppler frequency. Both sample movement and blood flow inside the tissue will change the reflected or scattered light frequency. However, changes caused by the movement of the sample and blood flow are quite different in nature. The blood flow in human beings is very complex due to the pulsatile nature, the parabolic velocity distribution and the randomness of the flow direction etc. In addition, cone-geometrical focusing beam, Brownian motion and the speckle will cause the broadening of the spectrum of the Doppler frequency shift and increase the variance. Sample movement is, however, much simpler. Most of bulk-motion correction methods assume that the sample movement will induce a constant phase term for each A line and by subtracting this phase term, the bulk motion induced artifacts can be corrected [15-17,23-25]. This assumption has proved to be valid and has been adopted by several groups to demonstrate great success [15-17,23-25]. Here, we also take this assumption.

If in Eqs. (2),(4)-(8) is the complex signal when there is no sample movement and is the bulk phase induced by sample movement at the A line, when there is sample movement, we can rewrite Eq. (2) as

where and is the Doppler frequency shift induced by the axial direction sample movement. So when there is axial sample movement, the average Doppler frequency shift detected is the linear sum of the Doppler frequency shift induced by sample movement with the Doppler frequency shift without sample movement.

As for the Doppler variance, we can rewrite Eq. (4) as following when there is sample movement,

We can find that the Doppler variance is actually only a function of the amplitude of the signal and not related to the phase term of the signal. In the case where lateral averaging is used as in Eq. (6), we can rewrite it as following when there is sample movement,

where and is proportional to the velocity of the axial direction sample movement. The Eq. (11) can also be written as

where and is proportional to the velocity of the axial direction sample movement. So, when there is sample movement and Eq. (6) is used to calculate the Doppler variance, the value of the variance is a function of , which is proportional to the sample-movement-induced velocity difference between the 1st A line and the other A lines in the averaging area. For most in vivo applications, the velocity induced by sample movement changes gradually and this velocity difference in the small averaging area (usually 4~16 A lines) is small. When both lateral and depth averaging algorithms are used as in Eq. (8), the impact of sample movement on the Doppler variance is reduced further (by N times) because is assumed to be constant in a single A line [15-17,23-25]. In this manuscript, we use Eq. (7) and Eq. (8) for the calculation of color Doppler and Doppler variance and , is used to balance the SNR, resolution and computation time.

4. Bulk-motion insensitive Doppler variance imaging

Details of the system used in this study can be found in references [17,26]. Briefly, the spectrometer-based Fourier domain OCT uses a Super luminescent Diode (SLD) light source which has a central wavelength of 890 nm and FWHM bandwidth of 150 nm. A modified scanning head from a commercial Zeiss Stratus OCT was used. The optical power on the human eye was set at 700 μW. The CCD integration time was set at 50 μs. The system sensitivity was measured to be about 100 dB at around zero imaging depth. The 6-dB sensitivity roll-off distance was found to be at an imaging depth of 1.6 mm. The imaging process includes background signal subtraction, linear interpolation to convert data from the linear wavelength space to the linear wavenumber space, and numerical dispersion. The axial resolution was measured to be 3.5 μm. The phase stability of the system is measured to be 0.48 milliradians.

As shown in Eq. (9), the detected mean Doppler frequency is the linear sum of the mean Doppler frequency shift of the blood flow with the Doppler frequency shift induced sample movement when there is sample movement. In order to accurately get the Doppler frequency induced by the blood vessel only, the sample induced bulk-motion must be removed. Up to now, several methods have been proposed to determine the bulk phase and most of the methods use a histogram-based method [15-17,23-25]. Usually, the median phase or mean phase are determined from the histogram of adjacent A line phase differences and then subtracted. Recently, averaged shift histogram based nonparametric density estimators have shown good performance and have become popular in extracting the bulk phase [15-17,24,25]. These histogram based methods have shown success in different situations; however these methods increase the computational complexity and time.

Figure 1 . shows the OCT images of the human retinal region. The images shown have a scanning range of 5 mm and include 2048 A lines. Figure 1 (a) shows the OCT structural image. Figures 1(b) and (c) show color Doppler images with and without bulk-phase correction. Figures 1(d) and (e) show the Doppler variance images with and without bulk-phase removal. The bulk-phase was removed by the histogram method introduced in reference 17. For the color Doppler and Doppler variance images, an intensity threshold that is 15dB more than the OCT intensity noise floor has been used to eliminate the low scattering but high phase noise region. For the color Doppler images, a threshold corresponding to a phase difference of 0.07 radians is also used to reduce the background noise but also consequently reduced the minimal detectable velocity with the color Doppler method to 99 μm/s. For the Doppler variance images, a threshold corresponding to a normalized Doppler variance of 0.05 is used to further reduce the background noise. All the color Doppler and Doppler variance images are processed based on Eq. (7) and Eq. (8).

Fig. 1

(a) OCT structure image; (b) Color Doppler OCT image without bulk motion correction; (c). Color Doppler OCT image with bulk motion correction; (e) Doppler variance OCT image without bulk motion correction; (f) Doppler variance OCT image with bulk motion correction. Scale bar: 1mm.

From Fig. 1(b), it can be seen that color Doppler is affected by bulk-motion. The bulk-motion increases the background signal of the color Doppler image and blood vessels cannot be indentified if the bulk-motion is too strong as shown in the region inside the red circle in Fig. 1 (b). It can also be found that the bulk-phase induced by bulk-motion is constant along the axial direction, which is the basis of the current bulk-motion algorithm. This also proves that the Doppler frequency induced bulk-motion can be considered as a constant variable in a single A line. Sample movement induced motion artifacts must be corrected before applying the color Doppler algorithm. The bulk-motion-corrected image as demonstrated in Fig. 1 (c) shows great improvement, and the color Doppler demonstrates much clearer blood vessels. However, in the region where large blood vessels exist, the bulk-motion-corrected image in Fig. 1 (c) shows correction artifacts as indicated by the red arrows. Although this correction artifact is not important in most imaging areas, in regions with large blood vessels such as the optic disk, the artifacts may cause erroneous blood vessel locations in the color Doppler image. An improved phase-resolved algorithm has been proposed to correct this artifact [26]. Although effective, the improved algorithm in reference 26 increases the calculation time.

In the Doppler variance image without bulk motion correction as shown in Fig. 1(d), the blood vessels can be seen even in the high speed sample movement region as shown in the red circle. As shown in the previous section, the Doppler variance obtained with the averaged autocorrelation algorithm is affected by the sample-movement-induced velocity difference within the averaging area. We can find from Fig. 1(b) that velocity of sample movement changes gradually (please note the green to red color change in the high speed sample movement region is due to phase wrapping). The sample-movement-induced velocity change is small in the small averaging area. Because depth averaging is also used, the bulk motion effect can be neglected here. The Doppler variance image with bulk motion correction is also shown in Fig. 1(e) in order to verify the results. The two images show great similarity even in the high bulk motion region as in the circle regions. However, just like in the color Doppler images, histogram-based bulk-motion-correction also introduces artifacts in the final Doppler variance image as indicated by the red arrow in Fig. 1(e). The Doppler variance image [Fig. 2(d) ] without bulk-motion-correction is artifact free.

Fig. 2

Speed-up curve for a workstation with dual quad-core CPUs of Intel Xeon processors X5550.

5. Real-time processing and displaying

Processing the Doppler OCT required more computation time than traditional structure OCT. Usually, data is acquired into memory and saved onto a hard disk, then processed afterward, or only a portion of the acquired data is processed and displayed in real time. For ophthalmological applications, real time processing and displaying of all the images is important. The commonly used histogram-based bulk-motion-removal method statistically analyzes the phase distribution along each A-line and extracts the bulk-motion phase. This method is time consuming and prevents real-time processing and displaying of Doppler OCT images. Doppler variance images are bulk-motion-correction free and artifact free in most applications as shown in the previous section. We have demonstrated a method to speed up the real-time processing time of polarization sensitive OCT by approximately 3 times with minimum modification of the existing program [27]. In this paper, we used a Dell workstation with dual Intel Xeon Processors X5550. Each Intel Xeon Processor X5550 has 4 cores and with a current dual CPU configuration, the system can support 8 threads running in parallel. Figure 2 shows the speedup curve of the current configuration. More than 7.5 times speedup was achieved compared with a single thread processing program. For the OCT system used in this study, the integration time for the CCD was set at 50 . The maximum frame rate was 9 frames per second when there were 2048 A-lines per frame. With a single thread involved in the data processing, the total processing time for OCT and Doppler variance images (without bulk-motion correction) was around 0.295 seconds. The Doppler variance images were processed with Eq. (8) with , . For the current configuration, 7 cores were dedicated to image processing and a speedup of 6.45 was achieved. For an image of 2048 A-lines, the processing time for OCT and Doppler variance images will take around 0.046 seconds. So the current system can process and display 20 frames per second in real time.

6. Results

In this paper, we were interested in the blood vessels in the choroidal layer. The retinal pigment epithelium (RPE) layer was first indentified based on the maximum OCT signal intensity along the axial direction [28]. The layers below the RPE layer were taken as choroidal layers. The RPE layer was identified based on the OCT structure images. The detected RPE layer location was applied to the Doppler variance images and color Doppler images. The Doppler variance images were processed without bulk-motion removal. The color Doppler images were post-processed with bulk-motion removal. The projection OCT structure, Doppler variance and color Doppler images were obtained by summing up all the layers below the RPE layer. The results are shown in Fig. 3 . The images are processed based on 120 frames of B-Scan images with 2048 A-line per frame and are acquired in around 12.3 seconds. The processing software is a homemade C++ program running on top of Microsoft Windows. Post processing is performed on a different computer from the one used for OCT system data acquisition. It took 3 minutes and 30 seconds to finish the whole processing and obtain the projection images. Different blood vessel network information can be obtained from these images. OCT structure images are based on the magnitude of backscattered or back-reflected light. Blood vessels usually manifest with a low magnitude in OCT structure projection images. Only large diameter blood vessels that are highly absorptive can be seen from OCT structure projection images. Color Doppler projection images [Fig. 3 (c)] can provide much more information about the blood vessel, such as blood flow direction (pink and blue colors represent different directions), blood flow speed, etc. However, due to the pulsatile nature of human blood flow, the blood flow direction and blood flow speed are not constant inside a single blood vessel. The incident angle dependent nature of Doppler OCT makes the situation more complex. The information provided by color Doppler OCT is compromised. Doppler variance projection image [Fig. 3 (b)] provides a fine microvascular network image which is more continuous and results in a better mapping of the vascular network. Because Doppler variance projection images are obtained without bulk-motion removal, the processing time is reduced and the processing procedure is much easier.

Fig. 3

(a)OCT projection, (b) Doppler variance projection and (c) color Doppler projection images of choroidal vessels. Scale bar: 1mm.

7. Conclusions

We demonstrated high resolution images of choroidal capillary vasculature network with phased-resolved color Doppler and Doppler variance. The effect of sample movement on the color Doppler and Doppler variance images was analyzed. It was found that the Doppler variance method is not affected by the sample movement. However, when the autocorrelation algorithm with averaging is used, the Doppler variance is affected by the sample-movement induced velocity changes among the different A lines in the averaging area and the effects of sample movement is reduced if the averaging is performed in both lateral and depth directions. When the sample-movement-induced velocity changes gradually, the Doppler variance images may be processed without bulk-motion-correction. Real-time processing and displaying the OCT and Doppler variance images with a parallel processing technique was also demonstrated.