1. Introduction

The Doppler effect was discovered by Christian Johann Doppler in the mid 1800s and has been widely used in many applications including medical imaging, particularly in the field of ultrasound imaging [1–3]. Analogous to ultrasonography, Optical Coherence Tomography (OCT) is a non-invasive optical imaging technology that is capable of depth sectioning of biological tissue at micrometer scale resolution [4]. Doppler Optical Coherence Tomography (DOCT) is an extension of conventional OCT that is capable not only of structural mapping in biological samples but also real time monitoring of flow characteristics such as location, direction, speed, and profile associated with the samples being imaged [5–7]. Various OCT-based flow detection techniques have been actively developed over the past ten years, such as the short time Fourier transform technique [8], the phase-resolved technique [9], optical microangiography (OMAG) [10], and time-frequency analysis called joint Spectral and Time domain OCT (STdOCT) [11]. Nevertheless, because of its simplicity in acquisition as well as in processing, Phase-Resolved DOCT (PR-DOCT) is the most widely used technique. PR-DOCT measures the amount of phase shift between two consecutive axial lines acquired at the same location in the sample and hence relies on the accuracy and stability of the detection of the phase difference. Knowing this phase difference Δϕ, the axial flow velocity can be determined as

 

Fig. 1 (a) Illustration of the incident beam angle in relation to the flow direction; (b) Fiber-based swept-source OCT built on Mach-Zehnder interferometer; PD = photo-detector, Cir = circulator, and PC = polarization controller.

Download Full Size | PDF

The phase difference Δϕ is normally computed through the use of the inverse tangential function and hence is subject to the π phase ambiguity. Therefore, without phase unwrapping, the maximum detectable phase shift is π, and the maximum detectable axial velocity is given by

The theoretical limit for the minimum detectable velocity, on the other hand, is determined by the phase stability of a PR-DOCT system. Different PR-DOCT systems exhibit different Doppler phase stabilities as characterized by a phase error in the absence of axial motion. This phase error determines the minimum detectable phase change and hence limits the velocity sensitivity of each particular PR-DOCT system. The phase error Δϕ_err can be statistically quantified as the deviation from the mean of the phase difference measured from a stationary sample. Consequently, the minimum detectable axial velocity (V_axial,min) is given by

The range from the minimum to the maximum detectable velocities, as defined above, determines the axial Velocity Dynamic Range (VDR) of each PR-DOCT system. While in principle, the upper limit can be extended by increasing the acquisition rate (i.e. 1/T), the later would cause the lower limit to also increase, which may be undesirable as this may lead to the invisibility of a slow flow [13]. Similarly, the lower limit may be extended by decreasing the acquisition rate but the upper limit would be also decreased, which would limit flow velocity measurement applications [7,14,15]. On the other hand, the lower limit may independently be lowered through the improvement of the stability of the system, which may be achieved in hardware but also in software by tailoring the algorithm used to calculate the phase difference [16]. Both the upper and lower limits of the VDR depend on the flow angle (θ) when the absolute flow is concerned. In practice, when monitoring an in vivo flow in a biological sample, the flow orientation may vary from zero to 90 degree relative to the incident beam, which leads to a wide dynamic range of axial flow velocity. It is to be noted that the invisibility of a slow flow will be more severe when the flow angle approaches 90 degree (i.e. flow parallels to sample surface), which produces extremely slow axial flow.

Recently, several alternative scanning protocols for different techniques of OCT-based flow detection have been developed to improve the sensitivity to slow flow of capillary blood vessels. The successful implementation of a method to increase flow sensitivity to an extremely slow flow was first demonstrated by Vakoc et al. [17]. Designing for three-dimensional (3D) flow segmentation, multiple axial scans (z-dimension) were acquired while performing two-dimensional (2D) scanning of the sample beam, consisting of fast scanning along the x-dimension and slow scanning along the y-dimension. To achieve high sensitivity to the slow flow, the phase differences were computed along the slow scanning direction. The flow segmentation was performed by computing the amplitude-weighted circular variance of the measured phase differences. The proposed technique required high oversampling (i.e. the ratio between the sample beam width and the lateral sampling interval) along the slow axis. The capability of the technique for in vivo 3D imaging of capillary vessel network over a wide field of view in mouse brain was demonstrated.

About the same time, Grulkowski et al. proposed several lateral scan protocols for achieving both fast-flow and slow-flow detections in a single frame acquisition implemented in their own technology called STdOCT [18]. STdOCT acquires multiple spectral interference signals over time at approximately the same lateral beam position and then extracts depth-resolved flow information by analyzing the amplitude shifted after 2D Fourier transformation of the acquired 2D spectral interference pattern. As a result, the proposed scanning pattern was designed for 2D oversampling as required by the processing algorithm of STdOCT. The technique required a tremendous amount of acquired spectra for each frame acquisition, which however was compensated by an ultrahigh speed spectrometer utilizing a high speed Complementary Metal-Oxide Semiconductor (CMOS) camera (i.e. 200 kHz line-rate). The detection of both fast-flow and slow-flow information in a single frame acquisition was proven to be useful for enhancing the visualization of the retinal capillary network in a flow segmentation application.

Most recently, an idea analogous to that presented in [17] was implemented in OMAG, an emerging technology that is capable of in vivo high resolution optical angiography [19,20]. Conventional OMAG introduces a constant frequency modulation along the fast scanning direction (i.e. x-direction) that causes the back-scattering signals from the moving and static scatterers to be separated after applying the OMAG algorithm [21]. The velocity sensitivity of OMAG is determined by the time interval between consecutive axial scans and hence the high-speed imaging requirement for 3D angiograms limits the sensitivity to slow flow in conventional OMAG. To overcome this limitation, a scanning protocol was modified and the OMAG algorithm was applied along the slow scanning axis (i.e. y-direction), which dramatically improved the slow flow detection of OMAG [19]. The 3D flow segmentation provided by OMAG is promising for the visualization of the 3D structure of the microcirculation of a blood vessel network.

In contrary to STdOCT and OMAG, which extract flow information based-on amplitude information, PR-DOCT extracts flow information directly from the phase shift between two axial scans acquired at different times. Nevertheless, high-speed PR-DOCT is also subject to limitation in the sensitivity to slow flow. The velocity sensitivity of PR-OCT is not only governed by the time interval between the two signals but also the phase stability of the system over that time period as can be observed from Eq. (3). As a result, the method of measuring the speed of the flow along the slow scanning axis using PR-DOCT is challenging given the expected low phase stability between the acquisition frames. Furthermore, the method of flow analysis along the slow scanning axis requires 3D scanning together with oversampling of the acquisition frame, adding complexity in the scanning protocol and involving large amounts of data in both the acquisition and processing stages. It is possible to apply the scanning protocol proposed in [18] to PR-DOCT but the increase in phase noise in the slow scanning axis is a potential problem. The degradation in phase stability may overthrow the benefit of having a longer time interval for slow flow detection. Moreover, PR-DOCT is theoretically capable of extracting flow information directly from only two axial scans, which is promising for real time flow monitoring. Nevertheless, the raster scanning protocol, where the sample beam continuously moves during a frame acquisition period, is known to be subject to severe degradation in Doppler phase stability unless accompanied by high oversampling [22,23]. The high oversampling requirement of the method involves a large amount of acquired spectra and computational times, which can be challenging for an application that demands real time processing and display.

In this paper, we present an alternative acquisition scheme along with a processing technique for a swept-source-based PR-DOCT that effectively extends the sensitivity to the slow flow and hence the detectable VDR of a given PR-DOCT system. The technique was designed to aim for real time acquisition, processing, and simultaneous display of multiple Doppler images having different detectable VDRs. Specifically, we present a technique of multi-scale measurement of flow velocity, where a sequence of spectra (i.e. up to 8 forward spectra) is acquired at each lateral position of the sample beam along the fast scanning axis and the Doppler phase shift is computed between each pair of spectra having different time intervals between them, yielding simultaneously different ranges of detectable flow velocity. The implementation of the technique for simultaneous acquisition and display of two Doppler maps having two different detectable VDRs (i.e. fast-flow and slow-flow Doppler maps) is demonstrated. To generate the two Doppler maps, a minimum requirement for the number of processed spectra per lateral position is only three and they do not need to be consecutive. Therefore, the processing involves a small number of axial scans and allows real time simultaneous display of the fast-flow and slow-flow maps. Finally, we demonstrate examples of real time Doppler imaging of an African frog tadpole using the technique. The overall imaging speed was limited in the case of our system by data read-out speed but was optimized to achieve real time in vivo imaging speed. The technique is particularly useful for phase-resolved Doppler detection utilizing high-speed Frequency Domain OCT (FD-OCT) where the acquisition time is extremely short and hence the ability to visualize the slow axial flow is often limited.

2. Method for multi VDR DOCT

The system was implemented with a swept-source-based FD-OCT built on a fiber-based Mach-Zehnder interferometer (MZI) as shown in Fig. 1(b). The FDML laser (Micron Optics) has a central wavelength of 1330 nm with a sweeping range of ~143 nm and an average output power of 12 mW. The output from the light source was coupled to the fiber and then split by an 80/20, 1x2 fiber coupler. The portion of the beam with 80% power was delivered to the sample arm. The lateral scanning in the sample arm was performed by using a galvanometer beam steering (VM500, GSI Lumonics). The sample beam was then focused into the sample through a 20 mm focal length plano-convex spherical lens. The 20% portion of the beam was delivered to the reference arm, in which the Fourier domain optical delay line was implemented in order to compensate for dispersion mismatch [24]. The two beams were coupled back to the fiber passing through fiber circulators and then combined at the 50/50, 2x2 fiber coupler. The time-encoded spectral interference signal was detected by using a balanced photo-receiver and then digitized on one channel of a two-channel, high-speed, 12-bit-resolution analog-to-digital converter operating at 200 Msamples/s (NI PCI 5124, National Instrument). The detected signal was recalibrated to the linear frequency-space prior to Fourier transformation to the depth profile. In this system, the recalibration process was done using a time-frequency relation measured by an unbalanced MZI denoted by a dash box in Fig. 1(b). Simultaneously with the detection of the main interference signal, the calibration signal was detected by a second photo-receiver and then digitized on another channel of the two-channel acquisition device. The calibration curve was generated from the position of peaks, valleys, and zero-crossing of an interference signal measured by the additional MZI [25].

The FDML laser was operated at a frequency sweep rate of 44.6 kHz. To extract the Doppler phase shift with multi-VDR capability, a block acquisition technique was implemented, where the number of sampling points was set so that multiple spectra consisting of both forward and backward sweep signals were recorded at exactly the same lateral position as shown in Fig. 2 . Furthermore, to avoid complexity in data processing and hence accommodate for real time processing and display, the backward sweep signals were omitted and only forward sweep signals as indicated by yellow dash boxes in Fig. 2 were used in the Doppler phase shift calculation. Therefore, the time interval between two consecutive forward spectra was 22.4 μs corresponding to a theoretical maximum detectable axial velocity of about 11 mm/s.

 

Fig. 2 Example of acquired spectra consisting 5 sweep cycles.

Download Full Size | PDF

The block of acquired forward sweep signals was then chopped into N sub-sections containing one forward spectrum per section. Each chopped signal contained 2000 sampled points. After calibration to the linear frequency-space, the number of sampling points per spectrum was approximately 1000 points. A fast Fourier transform (FFT) was performed with zero padding to 2048 points to increase the sampling resolution in the depth profile. Based-on Kasai autocorrelation [26], a technique that is commonly used in phase-resolved Doppler imaging [12,16,27], the phase shift was then calculated as a function of two variables using

Note that because Δϕ (z;p) is linear with p (i.e. Δϕ (z,p)/p is a constant for a given flow velocity), the detected velocity V _axial is only a function of z. Combining this calculation with the block acquisition technique, the VDR can be varied by changing the p parameter, which is equivalent to changing the time interval between two axial lines used to calculate the Doppler phase shift as illustrated in Fig. 2. It should be noted that, at the lower limit of the VDR as defined by Eq. (3), the detectable minimum phase shift is governed by the phase stability of the system that may vary as a function of p as will be shown in section 3. Therefore, to fully utilize the advantage of the multi-scale approach for slow flow sensitivity, the variation of the phase stability over p needs to be minimized so that the minimum detectable velocity decreases as p increases. At the upper limit of the VDR, the maximum detectable velocity is set by the π phase ambiguity that remains constant over the variation of p. By plugging T = pT ₀ in Eq. (2), the maximum detectable velocity therefore decreases as p increases. As a result, decreasing the time interval T alone will not extend the detectable VDR of the system. The technique presented here, which is capable of multi-scale measurement (i.e. simultaneous measurement of multiple VDRs), is required to effectively extend the overall VDR. The block acquisition technique was designed to address these two issues.

3. Phase stability of the multi-VDR DOCT

We demonstrate the experimental implementation of this formalism with the swept-source based PR-DOCT system capable of varying the VDR in real time, which is useful for the detection of the flow activity that exhibits a wide dynamic range of velocity. In this particular example, the number of sampling points was set so that each single block of acquired data contained eight forward sweep spectra (N = 8). The Doppler phase shift was then quantified by using Eq. (4) with M = 2 and p as a free-parameter, i.e. p = 1, 2, 3… 7, that could be changed in real time. The Doppler phase error corresponding to each value of p was quantified from the same set of acquired spectra by measuring the Doppler phase shift of a stationary mirror over time [13,28]. Examples of computed Doppler and intensity images of the stationary mirror are shown in Fig. 3(a) and 3(b), respectively.

 

Fig. 3 Illustration of the Doppler phase stability measurement; (a) Doppler image; (b) Intensity image; (c) zoom-in of the intensity image in (b); (d) axial profile of the intensity image in (c).

Download Full Size | PDF

Each image consisted of 100 axial lines (A-lines) corresponding with 100 lateral positions of the sample beam with a 10 μm sampling interval covering a lateral scanning range of 1 mm across the mirror surface. For each computed Doppler image, the measured phase shifts within the region around the mirror surface, as indicated in the red dash box in Fig. 3(c) (i.e. 4 pixels axially and 100 pixels laterally), were extracted and averaged. The number of pixels along depth was chosen to include only the region approximately within the Full Width at Half Maximum (FWHM) of the intensity axial profile around the signal peak as shown in Fig. 3(d), and the number of pixels along the lateral dimension covered the 1 mm lateral scanning range. The measurement was repeated 500 times (i.e. 500 Doppler images) and the histogram distribution was calculated as shown in Fig. 4 for the case of p = 1 and p = 6 as examples.

 

Fig. 4 An example of histogram distribution of the measured phase shift along with the corresponding Gaussian fitted curve (red dash line) for the case of (a) p = 1 and (b) p = 6; The horizontal axis is the phase shift error in mrad.

Download Full Size | PDF

It should be noted that, in this measurement, an intensity signal-to-noise ratio (SNR) was sufficiently high so that the effect of the background phase noise was negligible and hence the measured phase error represented the phase stability of the system as governed by, for example, the swept source, the interferometer, the scanning mechanics, and the signal processing. The SNR was computed as SNR = (S-μ_noise)²/σ ² _noise , where S was the measured signal peak amplitude at the position of the mirror surface averaged across the full lateral dimension of the acquired Doppler image, and μ_noise and σ_noise were the mean and standard deviation of the noise floor measured within the region around the signal peak while the sample beam was blocked, respectively [13]. The SNR in dB was defined by SNR_dB = 10log(SNR). Following the analysis in [22], the theoretical minimum detectable phase difference between two signals was determined as σ_{Δ ϕ} (rad) = 1/(SNR)^-1/2. An average SNR_dB at the signal peak was measured to be about 70 dB in this experiment, yielding a theoretical minimum detectable phase shift of about 0.3 mrad. Given that we will show thereafter that the measured phase error Δϕ _err is in the order of 10 to 20 mrad, the SNR_dB of 70 dB, which yields a theoretical phase error of 0.3 mrad, does not limit the phase stability measurement.

The FWHM of the histogram distribution was calculated representing the Doppler phase stability at each p value, as summarized in Table 1 . It should be pointed out that the phase error was computed from only two axial lines without averaging (M = 2) at each p value, and therefore represented extreme limit of the system phase stability within the time duration pT₀. Increasing the number of axial lines for each Doppler phase shift calculation, i.e. M > 2, will further improve the phase stability at the cost of reduction in imaging speed [6,16]. The phase stability tended to decrease as p was increased as shown in Table 1. However, since the time duration between the two axial lines used to determine the phase shift also increased, the minimum detectable velocity decreased as a function of p, and hence the method extended the total detectable velocity range of the system. From Table 1, the case of p = 1 represents the conventional PR-DOCT method, in which the detectable axial velocity is ranging from 46.5 µm/s to 10.9 mm/s, yielding a ratio between the maximum and minimum detectable velocity of about 230. Using the multi-VDR technique with N = 8, the detectable axial velocity is ranging from 11.3 µm/s to 10.9 mm/s that yields a ratio between the maximum and minimum detectable velocity of about 970, which is about four times of that achieved by the conventional method.

Table 1. Quantification of the detectable velocity range according to each value of p

View Table

4. In vivo imaging with multi-VDR DOCT

The application of the multi-VDR algorithm for real time in vivo imaging of biological sample is now detailed. The imaging scheme was designed in the way that two Doppler images were simultaneously acquired, processed, and displayed. One was determined at the shortest time interval between two spectra corresponding with the highest sweep rate of the light source that was capable of monitoring an axial flow speed of up to 11 mm/s in the sample. The other was determined at a longer time interval between two spectra (p > 2) that improved the sensitivity to the slow flow activity. The time duration T used to compute the second Doppler map was designed to be adjustable in real time by varying the value of p. All spectra needed to determine both Doppler images were acquired in a single shot at the same lateral position of the incident beam. We have shown in section 3 that the phase stability of the system was sufficiently high over different values of the p parameter. Therefore, in this experiment, the Doppler phase shift was extracted from only two spectra for each depth profile of the Doppler phase shift. For two Doppler maps (fast-flow and slow-flow) of size N _z x N _x pixels (depth x lateral), the number of processed spectra was 3N _x spectra. With such small amount of processed spectra, the system was capable of real time simultaneous display of the two Doppler maps along with the intensity map at a frame rate of about 3 to 4 frames per second. Each imaging frame consisted of 500 x 200 pixels (axial x lateral). The acquisition of the two Doppler images was considered as simultaneous because the acquisition time interval between them was less than 200 μs. The technique provides Doppler detection with extended velocity dynamic range as compared with a conventional PR-DOCT.

The blood flow activity within the heart of an African frog tadpole was chosen for demonstration of the multi-VDR DOCT technique. Structural images along with Doppler maps computed with p = 1 (VDR1) and p = 4 (VDR2) are shown in Fig. 5(a) –5(l). The images were taken at a certain cut plane around the location of the tadpole’s heart. The VDR1 alone is equivalent to the Doppler image that can be detected by the conventional PR-DOCT. It can be observed that certain speed of slow flow fell below the minimum detectable threshold of the VDR1 as demonstrated in Fig. 5(b). This slow flow however was picked up by the VDR2 as clearly observed in Fig. 5(c). Nevertheless, at maximum flow speed, the VDR1 provided un-ambiguity detection of the flow speed while the VDR2 was subject to the π phase ambiguity as shown in Fig. 5(k) and 5(l), respectively. Having both VDRs simultaneously therefore provides an extended range of detectable flow speed for Doppler imaging in a biological sample.

 

Fig. 5 (Color online) Doppler images of an African frog tadpole that were simultaneously acquired at the same location at different times where (a, d, g, and j) are intensity images, (b, e, h, and k) and (c, f, i, and l) are Doppler images corresponding to p = 1 and p = 4, respectively.

Download Full Size | PDF

Furthermore, we demonstrate the capability of varying the detectable velocity range in real time while monitoring an in vivo flow activity of an African frog tadpole using the multi-VDR technique (Media 1). A representative frame of the video recording of multi-VDR DOCT in action is shown in Fig. 6 . In this application, eight consecutive spectra were acquired simultaneously but only two spectra were processed at each lateral position. In reference to the Graphic User Interface (GUI) shown in Fig. 6, the time interval between two spectra used to compute the Doppler phase shift was changed by changing the value of p from 1 to 7 and the value of p + 1 is represent on the knob scale from left to right of the VDR control knob as demonstrated in Media 1. An image located on the upper-right corner is a conventional OCT intensity map representing the sample structure. On the lower-right corner is a Doppler image revealing flow information within the heart of the tadpole. The Doppler map shown on the GUI was captured at p = 3, which provides a mapping of bidirectional axial flow velocity ranging from −3.6 mm/s to + 3.6 mm/s with axial velocity sensitivity of 18.5 μm/s (see Table 1).

 

Fig. 6 (Color online) A representative frame of a video recording of flow activity within a tadpole heart acquired by multi-VDR demonstrating the capability of varying the detectable VDR of the Doppler map in real time (Media 1).

Download Full Size | PDF

In practice, even in the same sample, the flow in different regions can exhibit significantly different values in velocity dynamic range caused by, for example, flow orientation, vessel size, and functionality. By increasing the p value in multi-VDR DOCT, even though the maximum detectable velocity is limited by the phase wrapping, it is capable of revealing more flow details because of the high sensitivity to slow flow as illustrated in Fig. 7 . The intensity map provides depth cross-sectional of the area of interest around the heart chamber of the tadpole as shown in Fig. 7(a). The Doppler image in Fig. 7(b) maps an axial flow velocity within a non-ambiguity velocity range from −10.9 mm/s to + 10.9 mm/s with the velocity sensitivity of about 47 μm/s. The positive and negative flow velocities represent flow in opposite directions (i.e. positive velocity corresponds to flow moving toward the observer). An axial flow speed of about −6 mm/s is observed at the location pointed by the white arrow in Fig. 7(b). Simultaneously, the Doppler image in Fig. 7(c) provides flow information within an axial flow velocity range from −2.7 mm/s to + 2.7 mm/s with the velocity sensitivity of about 16 μm/s. Additional flow information is clearly visible in this Doppler map as indicated by the white arrow in Fig. 7(c). The axial velocity within this slow flow region is estimated to be less than −1 mm/s and hence cannot be observed in Fig. 7(b). Nevertheless, the fast flow region (i.e. ~6 mm/s axial flow speed) as indicated in Fig. 7(b) is subject to phase-wrapping in Fig. 7(c) and hence in the slow-flow map, the fast flow velocity cannot be determined without phase unwrapping. Therefore, the technique of multi-scale VDR is useful for real time monitoring of the flow activity in a biological sample.

 

Fig. 7 (Color online) (a) Conventional intensity image, (b) fast-flow Doppler image, and (c) slow-flow Doppler image acquired at the same location of the tadpole heart region.

Download Full Size | PDF

The main factor that limits the current frame rate is the data read-out time (i.e. the time for transferring the acquired data from the acquisition device to computer memory), which was long compared to the actual time for data sampling. The sampling time of one block of data consisting of 8 forward spectra took only about 180 μs and hence 36 ms was required per acquired frame. The average time for data transfer was approximately 200 ms per frame. Nevertheless, there is still room for improving the frame rate. To increase the imaging speed in the future, a more efficient design of the hardware interface as well as the synchronization in the data acquisition will be investigated. In this developing stage, the system was solely implemented in a standard Labview programming. The improvement in overall acquisition time as well as in the processing time can be achieved through the implementation of the technique in low level software programming such as the ANSI C programming or direct hardware programming such as a field programmable gated array [29].

5. Conclusion

In conclusion, while the maximum detectable flow speed of PR-DOCT can be extended through the shortening of the acquisition time interval, the minimum detectable velocity will also be increased; this approach leads to the invisibility of the slow flow activity whose axial speed is below the minimum detectable threshold of the system. We presented a technique of multi-scale flow velocity measurement that simultaneously acquires, processes, and display multiple Doppler images having different ranges of detectable flow velocity. The technique effectively extends the overall detectable velocity dynamic range of PR-DOCT. The Doppler phase stability and the corresponding minimum detectable velocity at different time intervals (T_p) were quantified. A detectable axial velocity ranging from 11.3 µm/s to 10.9 mm/s was achieved. Finally, we demonstrated the implementation of the proposed technique for in vivo Doppler imaging of an African frog tadpole. The results demonstrated the extended sensitivity to the slow flow activity that could not be detected by the conventional method.