1. Introduction

Doppler optical coherence tomography (DOCT) is a powerful technique for determining the localized flow velocity in highly scattering media such as biological tissues [1]. Axial component of the velocity vector is directly detectable in a DOCT, but the other two components of the 3-D velocity vector have been ascertained through ingenious approaches. Dave and Milner [2] were the first to measure Doppler angle by detecting velocity at two different angles with orthogonal polarization modes, where two Wollaston prisms were adopted with a prerequisite of polarization-insensitive property of the sample. Pedersen et al [3] used a glass plate to create delay-encoded sub-beams with two incident angles for determining Doppler angle. Doppler spectrum broadening [4] was exploited to measure transversal velocity and hence the Doppler angle. However, the azimuth angle of the velocity vector was still undetermined. Full determination of velocity vector was realized by a four-channel quadrant detector in a free-space system [5], but such implementation is not applicable to a fiber-based system. A beam divider was proposed to quantify the 3-D velocity vector in a fiber-based spectral-domain DOCT system [6]. The beam divider encoded the probe beam with five different path length delays requiring too much depth range, which causes severe restriction on spectral-domain OCT system. Other groups [7–9] tried to get the velocity vector from the estimation of the flowing direction based on reconstructed 3-D orientation of a blood vessel. However, the velocity direction of different position inside the vessel is not strictly consistent with the vessel orientation.

In this paper, we present a transit-time based method in a fiber-based spectral-domain DOCT system to ascertain the azimuth angle of the velocity vector. The method is especially useful to determine the azimuth angle with high accuracy within 90° range under prior knowledge of azimuth angle known from the reconstructed 3-D orientation of the blood vessel. By inserting a slit plate with predetermined slit orientation into the sample arm, the probe beam is divided into three delay-encoded sub-beams with different beam shape. Three independent complex signals corresponding to different path length delays are obtained by Fourier transform of the detected spectral interferogram. The azimuth angle within 90° range is then determined by Doppler bandwidth measured from the three complex signals. 3-D velocity vector is thus quantified by further estimating of Doppler angle and flow velocity by combined Doppler shift and Doppler bandwidth measurements [10]. The available depth is 1/3 of the maximum imaging depth in the proposed method, in contrast to 1/5 of the maximum imaging depth in the previously reported method [6]. Therefore, the demand on imaging depth for spectral domain DOCT system is relaxed. Higher accuracy on determining velocity vector is hence envisioned due to the depth dependent sensitivity fall-off in spectral-domain OCT. Finally, feasibility of the proposed method is demonstrated by phantom experiments.

2. Theory

The method is based on the full exploitation of Doppler spectrum broadening. There are a number of sources, including flow velocity, turbulence, Brownian motion, as well as the probe beam geometry, contribute to the broadened Doppler spectrum [4]. Doppler spectrum broadening is mainly determined by the transit time of scattering particles passing through the focused probe beam [11]. The transit time in the focal area is

 

Fig. 1 Custom-designed slit plate for the quantification of 3-D velocity vector. Inset figures demonstrating the top view of the slit plate with specified parameters (a), and the layout of the sample arm in the spectral-domain DOCT system (b).

Download Full Size | PDF

To quantify the 3-D velocity vector of moving particles, it is required to obtain three quantities, which are flow velocity V, Doppler angle α and azimuth angle φ. As shown in Fig. 1(a), the slit plate is made of glass with air slit in the middle. The glass part of the plate has a thickness of t with refractive index of n_g. The slit with a width of d is oriented along x axis. The plate is introduced in the sample arm of a spectral-domain DOCT system between the collimator and galvo mirror shown in Fig. 1(b). Collimated beam with a diameter of D illuminates the plate and focused to the scattering particles by a focusing lens with the focal length of f. Backscattered light from the moving particles is collected by the focusing lens and passes through the plate again. In double passing of the plate, the light could go through the plate via air slit first and then return via air slit or glass. Alternatively, the light could go through the plate via glass first and then return via air slit or glass. Therefore, the slit plate makes four different beam paths: air to air (AA), air to glass (AG), glass to air (GA) and glass to glass (GG). Since AG and GA have the same optical path delay, there are three independent path length delays, which can be denoted by 0, t(n_g-1) and 2t(n_g-1). In the following text, the double-subscript of the transit length w and the Doppler bandwidth B indicates the corresponding beam path.

The maximum imaging depth in optical length of the spectrometer based spectral-domain DOCT system is determined by Mλ ²/4Δλ, where M is the number of pixels of the CCD, λ is the central wavelength of the light source, and Δλ is the full spectral range of the spectrometer. In the previously reported method [6], neighbouring sub-images of the sample with an optical depth larger than Mλ ²/20Δλ will overlap, while in the proposed method this threshold is improved to be Mλ ²/12Δλ. Therefore, the demand on imaging depth for the spectral-domain DOCT system is relaxed.

Figure 2(a) and (b) shows the simulation results of normalized intensity patterns in the focal plane of the focusing lens corresponding to collimated illumination via air slit (d=1 mm) and glass, respectively. x and y are horizontal and vertical coordinates of the focal plane. To analyze the transit time, the beam size is estimated from the contour of the central maximum of the intensity pattern. Under illumination via air slit demonstrated in Fig. 2(a), the width of the central maximum along x direction and y direction are 2λf /D and 2λf /d, respectively. The fitted contour of the central maximum is almost a rectangle. Therefore, the transit length w _AA that a particle passes through the central pattern with an azimuth angle φ can be expressed as

 

Fig. 2 Simulation results of normalized intensity patterns in the focal plane corresponding to collimated illumination via air slit (a) and glass (b), respectively. Contours of the central maximum corresponding to case (b) with different width of air slit (c).

Download Full Size | PDF

Under illumination via glass as demonstrated in Fig. 2(b), contour of the central maximum can be fitted by elliptic equation. The fitted contours with different width of d are plotted in Fig. 2(c). Obviously, the central maximum under illumination via glass approaches an Airy pattern as the air slit d decreases. The transit length w _GG that a particle passes through this central pattern with an azimuth angle φ is given by 2(x ₀ ²cos² φ+y ₀ ²sin² φ)^1/2, where x ₀ and y ₀ represent the semimajor axis and semiminor axis of the fitted ellipse, respectively. For instance, with system parameters of λ=0.84 μm, f=75 mm, D=5.2 mm and d=1 mm, the fitted elliptic contour of the central maximum is (x/15.77)²+(y/11.50)²=1. Thus we have x ₀=15.77 μm and y ₀=11.50 μm.

The azimuth angle within 90° range is ascertained by full exploitation of complex signals corresponding to three path length delays. In spectral-domain DOCT, a depth profile is obtained by calculating the Fourier transform of the spectral interferogram in the k-space. The depth profile can be recognized as a complex signal represented by $\tilde{\Gamma }(z)=A(z)\exp [i\Phi (z)]$, where z is the axial coordinate, A(z) and Φ(z) denote the amplitude and phase terms, respectively. A(z) is used for structural information reconstruction, and phase Φ(z) is used for Doppler information construction. Then $\tilde{\Gamma }(z)$ is separated into three complex signals corresponding to three path length delays. Assuming ${\tilde{\Gamma }}_{AA}(z)$ is the complex signal corresponding to sub-beam AA, ${\tilde{\Gamma }}_{AA}(z)$ is extracted from $\tilde{\Gamma }(z)$ with z in the range of (z₀-δ/2, z₀+δ/2), where δ=t(n_g-1)/2 and z₀ represents the axial position corresponding to the center of sample depth. ${\tilde{\Gamma }}_{AG+GA}(z+\delta )$ is the complex signal corresponding to sub-beams AG and GA, which adds a delay shift of δ in axial direction compared with ${\tilde{\Gamma }}_{AA}(z)$. Similarly, ${\tilde{\Gamma }}_{GG}(z+2\delta )$ is the signal corresponding to sub-beam GG, which adds a delay shift of 2δ in axial direction compared with ${\tilde{\Gamma }}_{AA}(z)$.

By coherently summation of the complex signals with three path length delays, we have

According to Eq. (2), when the background B ₀ is subtracted, the ratio of the Doppler bandwidth obtained from ${\tilde{\Gamma }}_{AA}(z)$ to that obtained from ${\tilde{\Gamma }}_{sum}(z)$ is given by

Similarly, the ratio of the Doppler bandwidth obtained from ${\tilde{\Gamma }}_{GG}(z+2\delta )$ to that obtained from ${\tilde{\Gamma }}_{sum}(z)$ is given by

Based on the relationship between Doppler bandwidth B and standard deviation σ, i.e. $B=4\sigma$, the Doppler bandwidth is obtained from the standard deviation, which could be calculated from complex signals of successive A-scans according to [4]:

After calculating Doppler bandwidth B _AA, B _GG and B _sum, azimuth angle φ within 90° range can be determined using Eq. (5) or (6). However, the determination of an unknown φ through Eq. (5) is under a prior knowledge of its range, which can be obtained by Eq. (6). The other two quantities required for the quantification of the velocity vector, Doppler angle α and flow velocity V, are estimated by combined Doppler shift and Doppler bandwidth measurements from the coherently summed complex signal ${\tilde{\Gamma }}_{sum}(z)$ based on previously reported method [10].

The relationship between Doppler shift f_d and flow velocity V is

Doppler shift f_d is calculated by phase difference of successive A-scans:

From Eq. (2) and Eq. (8), the Doppler angle α is calculated by

The flow chart of signal processing procedure is shown in Fig. 3 .

 

Fig. 3 Flow chart of signal processig procedure to quantify the velocity vector.

Download Full Size | PDF

Intuitively, the case via air slit shown in Fig. 2(a) is more sensitive to slit orientation in comparison with that via glass shown in Fig. 2(b), and higher accuracy on determination of azimuth angle might be envisioned using Eq. (5) than using Eq. (6). However, the SNR corresponding to ${\tilde{\Gamma }}_{AA}(z)$ is worse than ${\tilde{\Gamma }}_{GG}(z+2\delta )$ due to the small width of air slit d. Therefore, both Eq. (5) and (6) are tried for the evaluation of azimuth angle in the following experiments.

3. Experiments

To validate the feasibility of the transit-time based method for 3-D velocity vector quantification, the fiber-based spectral-domain OCT system previously described [13, 14] is functionally extended to a DOCT system besides inserting the slit plate into the sample arm as shown in Fig. 1(b). The light source is a superluminescent diode with a central wavelength of 840 nm and bandwidth of 45 nm. In the sample arm, the 5.2 mm diameter beam passes through the slit plate, scanned by a galvo mirror and focused by a 75 mm focal length lens. In the spectrometer, the interfered beam is collimated by a 60 mm focal length lens, dispersed by a volume phase holographic grating of 1200 lines/mm, focused by a 150 mm focal length lens, and detected by a line-CCD camera of 2048 pixels with 20 KHz A-scan rate. Using a 1.2 mW sample beam, the sensitivity measured near DC is 108 dB with an axial resolution of 7.5 µm in air. The dynamic range is 68 dB. There is a 16 dB falloff in sensitivity over 2 mm depth. The unambiguous velocity range lies on the maximum detectable Doppler shift, which is restricted by the A-scan rate. The determinable Doppler angle is from 0° to 180°. The determinable azimuth angle is within 90° range inside a quadrant.

0.25% aqueous solution of polystyrene beads (0.53 μm diameter) in a capillary tube with inner diameter of 150 μm is used as the sample. Flowing speed of the scattering beads within the tube is controlled by a syringe pump. The tube orientation is set at azimuth angle φ starting from x axis and Doppler angle α starting from z axis as depicted in Fig. 1. The air silt width d is set to be 1 mm or 0.8 mm for comparable signal intensities of three sub-images corresponding to the three path length delays. The glass part of the plate is made of BK7, which has a thickness of 1.7 mm with a refractive index of 1.51 at 830 nm wavelength. In our experiment, group velocity mismatch among the different paths introduced by the slit plate is not considered since the thickness of the glass part is relatively small.

Figure 4 demonstrates the cross-sectional Doppler images representing Doppler shift and Doppler bandwidth of the capillary tube with flowing solution of polystyrene beads inside. The three sub-images of the same sample corresponding to three path length delays are obtained simultaneously by the fiber-based spectral-domain DOCT system.

 

Fig. 4 Doppler shift image (left) and Doppler bandwidth image (right) of the capillary tube with the slit plate inserted in the sample arm corresponding to parameter setting of α=76°, φ=45°, V=4.72 mm/s, d=1 mm.

Download Full Size | PDF

To investigate the ability of the proposed method to determine azimuth angle, experiments corresponding to ten different preset azimuth angles within 90° range are conducted with slit width d of 1 mm and 0.8 mm. Under the condition that the orientation of capillary tube at a fixed Doppler angle of 76°, the volume velocity of flowing solution is varied from 1 μl/min to 5 μl/min at increment step of 1 μl/min, equivalent to an average velocity changing from 0.94 mm/s to 4.72 mm/s.

M-mode data acquisition is to perform several A-scans at a transversal position without transversal scanning [4]. In the experiment, the probe beam is focused to the center of the capillary tube. N A-scans are performed for averaging (N=1000). With A-scan rate at 20 KHz, the time needed for one M-mode data acquisition is 50 ms. The reason why we use M-mode in the experiment is to exclude other potential affecting factors which include: 1) deflexion of incident ray caused by geometric shape of capillary tube and the refraction index mismatch; 2) beam shape variation with the position of the scanning beam due to optical aberrations; 3) additional phase noise induced by the transversal scanning.

After calculating the Doppler bandwidth obtained from ${\tilde{\Gamma }}_{AA}(z)$, ${\tilde{\Gamma }}_{GG}(z+2\delta )$ and ${\tilde{\Gamma }}_{sum}(z)$ under different flow velocities, the corresponding background B ₀ can be deduced through linear fitting and extending to zero velocity. The estimated background B ₀ for cases of ${\tilde{\Gamma }}_{AA}(z)$, ${\tilde{\Gamma }}_{GG}(z+2\delta )$ and ${\tilde{\Gamma }}_{sum}(z)$ with d=1 mm are 128.8, 140.4 and 136.0 Hz, respectively, and with d=0.8 mm, those are 121.6, 139.2 and 133.2 Hz, respectively.

The azimuth angles estimated by Eq. (5) and Eq. (6) with preset flow velocity of 4.72 mm/s are shown in Fig. 5(a) and (b) , respectively. As shown in Fig. 5(a), the estimated azimuth angles are in good agreement with the preset ones for φ smaller than 80°. The accuracy for two values of d is comparable. Invalid estimation of azimuth angle occurs when it approaches 90° due to the big influence of Doppler bandwidth background, thus the calculated data at 80° and 90° are discarded. Compared with the results shown in Fig. 5(a), the estimated azimuth angles in Fig. 5(b) are in better agreement with the preset ones especially for φ near 90°. The data for φ around 45° is in perfect match with the preset ones, since w _GG as a function of φ is approximately linear in this situation. However, the data for φ near 0° or 90° is not accurate for the reason that w _GG is insensitive to the variation of φ in these two situations. As can be seen, the result for d=1 mm is in better agreement with the preset ones compared with d=0.8 mm. As d becomes larger, the fitted elliptic contour of the central maximum under illumination via glass becomes more elongated (see Fig. 2(c)). Therefore the transit length w _GG is more sensitive to the variation of φ. The data of 0° is discarded since the calculated cos² φ (see Eq. (6)) is invalid owing to added noise. The RMS errors of the estimated angles corresponding to Fig. 5(a) and (b) for two values of d are summarized in Table 1 . The difference between estimated azimuth angles based on Eq. (5) and preset ones is within the range of 0.29° to 9.37°, while based on Eq. (6) the difference is within the range of 0.19° to 12.90°. Higher accuracy of estimation on azimuth angle is achieved using Eq. (5) than using Eq. (6) if invalid estimations at 80° and 90° are excluded. On the other hand, method based on Eq. (6) can still obtain effective estimation of azimuth angle even it approaches 90° due to relatively high SNR.

 

Fig. 5 Estimated azimuth angle using B _AA/B _sum (a) and B _GG/B _sum (b) with different preset azimuth angles for two values of d corresponding to parameter setting of α=76°, V=4.72 mm/s.

Download Full Size | PDF

Table 1. RMS Errors for Estimation of Azimuth Angle

View Table

After azimuth angle φ is determined, the Doppler shift and Doppler bandwidth obtained from the coherently summed complex signal ${\tilde{\Gamma }}_{sum}(z)$ are used for the estimation of flow velocity V and Doppler angle α. The estimated results are shown in Fig. 6(a) and (b) , where the capillary tube is positioned at azimuth angle of 45°. It can be seen that the estimated data for d=1 mm and 0.8 mm are both in good agreement with the preset ones under different flow velocities. The estimated flow velocity vector is thus quantified and plotted in Fig. 7 with the preset flow velocity of 4.72 mm/s.

 

Fig. 6 Estimated flow velocity (a), Doppler angle (b) with different preset flow velocities for two values of d corresponding to parameter setting of α=76°, φ=45°.

Download Full Size | PDF

 

Fig. 7 Reconstructed flow velocity vector profile inside the capillary tube corresponding to parameter setting of α=76°, φ=45°, V=4.72 mm/s.

Download Full Size | PDF

In the text above, the Doppler bandwidths are computed for in focus situation. In the situation out of focus, the actual transit length w is enlarged due to bigger lateral spot size. Therefore, according to the relationship of Doppler bandwidth and transit length in Eq. (2), the Doppler bandwidth B decreases. In our method, the azimuth angle φ is obtained by calculating the ratio of B _AA and B _sum or B _GG and B _sum according to Eq. (5) and (6). Although the Doppler bandwidth decreases, the variation of B _AA/B _sum or B _GG/B _sum is usually negligible since the decreasing of B _AA, B _GG and B _sum is almost proportional, especially when the defocusing distance is small.

However, defocusing effect will result in an incorrect estimation of Doppler angle. The Doppler angle is calculated by the combination of Doppler bandwidth and Doppler shift as demonstrated by Eq. (10). Where w is the theoretical transit length, f_d is the Doppler shift which is not affected by defocusing. Therefore, the Doppler angle is underestimated.

In the sample arm of our system, the 5.2 mm diameter beam is focused by a 75 mm focal length lens, thus the depth of focus is 445 μm. The depth of focus is much larger than the optical depth of the experimental phantom. The probe beam is weakly focused, and in this case, the defocusing effect is controlled to the minimum.

4. Conclusion

In summary, we present a method to quantify a velocity vector in spectral-domain DOCT based on transit-time analysis of scattering particles passing through the focusing pattern formed by delay-encoded sub-beams. Doppler bandwidth measured from the complex signals corresponding to three different delays as well as the coherently summed complex signal are fully exploited for the determination of azimuth angle. From conducted experiments on azimuth angle estimation of a capillary tube at preset orientations, we compare the accuracy of two methods based on Doppler bandwidth ratio measured from different complex signals. Velocity vector of flowing solution inside a capillary tube quantified with discernible velocity profile is in good agreement with expectation. Under prior knowledge of azimuth angle known from the reconstructed 3-D orientation of the vessel, it is feasible to determine the azimuth angle with high accuracy within 90° range and thus obtain velocity profile within the vessel. Further study on the numerical dispersion compensation is necessary in the application of highly dispersive biomedical sample.