Improved phase-resolved optical Doppler tomography using the Kasai velocity estimator and histogram segmentation

Signiﬁcant improvements are reported in the measurable velocity range and tissue motion artefact rejection of a phase-resolved optical coherence tomography and optical Doppler tomography system. Phase information derived from an in-phase and quadrature demodulator is used to estimate the mean blood ﬂow velocity by the Kasai autocorrelation algorithm. A histogram-based velocity segmentation algorithm is used to determine block tissue movement and remove tissue motion artefacts that can be faster or slower in velocity than that of the microcirculation. The minimum detectable Doppler frequency is about 100 Hz, corresponding to a ﬂow velocity resolution of 30 l m/s with an axial-line scanning frequency of 8.05 kHz and a mean phase change measured over eight sequential scans; the maximum detectable Doppler frequency is (cid:1) 4 kHz (for bi-directional ﬂow) before phase wrap-around. (cid:1) 2002 Elsevier Science B.V. reserved.

Developed by several groups for imaging blood flow non-invasively, optical Doppler tomography (ODT) [1][2][3][4][5][6], or color-Doppler optical coherence tomography (CDOCT) [7][8][9][10][11], is an extension of optical coherence tomography (OCT) [12]. The exceptional spatial resolution of OCT allows visualization of structures close to the cellular level in vivo [13,14], while the high velocity sensitivity of ODT permits imaging of the microcirculation in blood vessels too small to be resolved by conventional Doppler ultrasound. Previously, a phaseresolved OCT/ODT system has been reported [3][4][5][6], in which the phase information derived from the Hilbert transformation of the interference sig- 15 July 2002 Optics Communications 208 (2002) 209-214 www.elsevier.com/locate/optcom nal is used to obtain the blood flow velocity distribution in tissue. Unlike other methods that use short-time fast Fourier transformation (STFFT [1,2,[7][8][9][10][11], this technique permits increased axial scanning speed while maintaining good velocity sensitivity. However, deriving high-quality phase information by Hilbert transformation demands fast analog-to-digital conversion and processing of large amounts of digital data, resulting in long computation times that can be the limiting factor in achieving higher frame rates. Here, we describe an improved system using hardware in-phase and quadrature (I&Q) demodulator with software implementation of the Kasai autocorrelation velocity estimator [15]. The significant reduction in computation time not only permits in vivo imaging of blood flow (demonstrated here in human skin), but also allows determination of the orientation of the blood vessel by repositioning the probe interactively. Hence, the Doppler angle can be measured to determine the absolute blood flow velocity.
A crucial issue in imaging microcirculation is that the blood flow velocities can be orders of magnitude lower than typical physiological or bulk tissue motion. Rejection of the tissue motion artefact thus becomes critical. Previously, crosscorrelation and edge-detection registration algorithms were used to remove such artefacts [11,16] in retinal imaging applications. However, these methods are very dependent on prominent structural features such as the bright reflections or edges that may not be available in many tissues. Here, we present a histogram segmentation algorithm using features in the calculated velocity distribution function, with the assumption that ODT is usually employed to image small blood vessels, as explained later.
A schematic of the ODT system is shown in Fig.  1. A broadband 1.3 lm light source (JDS, Ont., Canada) with a polarized output power of 5 mW and bandwidth of 62 nm is used. The coherence length of the source is measured as $13 lm in air, corresponding to an axial resolution of $10 lm in biological tissue (assuming tissue index of refraction $1.4). A rapid-scanning optical delay (RSOD) line is used in the reference arm, similar to other high-speed OCT systems [17,18]. The RSOD is aligned such that no phase modulation is gen-erated when the group delay is scanned sinusoidally at 8.05 kHz. The phase modulation is generated using an electro-optical phase modulator (JDS, ON, Canada). In order to accommodate the high RSOD scan frequency, the carrier frequency generated by the phase modulator is set at 4.3 MHz. An optical circulator (O-Net, Boston, MA) is used to increase signal-to-noise ratio (SNR) [19,20]. A balanced photodetector of 10 MHz bandwidth to measure the interference fringes, and an I&Q demodulator of 1.55 MHz bandwidth to obtain the base-band in-phase signal I, and quadrature signal Q, were constructed. A reference 4.3 MHz sinusoidal signal synchronized to the phase modulator is provided to the demodulator. The SNR of the system is optimized such that the photon shot noise from the reference arm dominates over electronic noise. The phase drift is approximately p rad/s after the RSOD scanner is warmed up. I and Q are digitized using a high-speed 12-bit analog-to-digital converter at 5 MHz (National Instrument, Austin, TX) into a personal computer with dual Pentium III 850 MHz processors. Both the inward and outward portions of the depth scan are digitized. The sample arm contains an infrared laser-focusing doublet lens, producing a beam waist diameter of about 10 lm. A linear translator (Newport, Irvine, CA) capable of 6 mm/s maximum velocity is used for lateral scanning. For structural imaging, eight A-scan (i.e., depth scan) lines are averaged to improve the SNR such that a 312 Â 600 pixel image is acquired in 0.6 s. For Doppler imaging, these eight A-scan (i.e., depth scan) lines are used in Eq. (2) (see below) to estimate the velocity distribution.
At any particular depth, the complex digitized signal can be expressed as C ¼ I þ jQ. The structural OCT image is then computed for each pixel as A ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi . The mean phase change over N depth scans computed by the Kasai autocorrelation [21] algorithm is where N is the number of averaged A-scans and i denotes the signal of the ith A-scan at a given depth. A four-quadrant arctan function is used to compute the phase change in the interval [Àp; p].
The Doppler frequency shift induced by a moving particle is then f D ¼ D/=2pT , where T is the time between successive A-scans, and the mean velocity is estimated by where k 0 is the source center wavelength, n is the tissue index of refraction, and h is the angle between the direction of movement and the optical axis, i.e., the Doppler angle. The estimated mean velocities of all pixels in an A-scan line are then used to compute a histogram. Assuming the blood vessels being measured have small diameters compared to the A-scan depth, the velocity of any block tissue motion will appear at the peak of this histogram. The velocity at which the histogram reaches a maximum is then subtracted in the interval [Àp; p] to reject the tissue motion artefacts. For this experiment, a 128-bin histogram was used to estimate the block tissue motion velocity. In regions of low SNR, random noise from the I and Q channels dominates and produces erroneous velocity estimates. A threshold (SNR > 2) is, therefore, used to suppress such effects to produce the final color-Doppler flow image [21]. Without considering speckle effect, the minimum resolvable velocity is related to the minimum detectable phase change, measured to be approximately 3°. Fig. 2 shows the overall signal flow chart of the processing algorithm.
The system was used to image subsurface microcirculation in the finger skin of a human volunteer in vivo. The high A-scan frequency at 8.05 kHz and the tissue motion artefact rejection algorithm allowed these measurements to be performed without physically immobilizing the tissue. A thin layer of index-matching oil was applied to the surface of the skin to increase light penetration into the tissue. A structural OCT image and the corresponding color blood flow image are shown in Fig. 3. When searching for blood vessels, the frame time is $1 s including signal acquisition ($0.6 s) and data processing ($0.4 s) for both the structural and flow images (two frames are generated using the inward and outward A-scans respectively). After a blood vessel of interest was found, the lateral scanning speed was decreased to reduce the velocity noise caused by speckle modulation and to produce a high quality color flow image. In the center of Fig. 3(b), a capillary of approximately 20 lm diameter is seen with a peak Doppler frequency shift of about 640 Hz. This capillary extends for some 500 lm when sequential frames are viewed. Using the 1 frame/s image update, the subject moved the finger interactively such that a large portion of the capillary was brought into alignment with the imaging plane. The Doppler angle was then measured as 78°. By Eq. (2) the peak flow velocity at the center of the vessel was 1.4 mm/s, in good agreement with pre- vious measurements [4,22,23]. The noise in the velocity image is approximately 100 Hz, corresponding to about 5°phase noise at the 8.05 kHz A-scan frequency. This is slightly higher than the 3°phase noise when the system is not scanning laterally, and is likely due to the speckle modulation caused by the lateral motion and/or vibration caused by the scanning motor.
Coherent demodulation has been demonstrated previously [7,9] in OCT and CDOCT systems, but they mainly employed STFFT-based algorithms to obtain flow information. This has high computational complexity and requires $10 s per frame for offline processing [7,8]. High-speed resonant scanners have been used in OCT systems [20] to achieve high frame rate, but the carrier frequency of the signal in these systems results from phase modulation in the RSOD [18], which varies sinusoidally and can deviate due to the angular amplitude fluctuation of the scanning mirror. To avoid these difficulties in our current setup, the electro-optical phase modulator allows centering the diffracted light beam on the RSOD scanning mirror, giving a more stable and constant carrier frequency (drift < 0:5 Hz) at 4.3 MHz. This relaxes the demand of filtering in the I&Q demodulator even for a very high A-scan rate ($8 kHz) and allows the application of the Kasai estimator. The higher A-scan rate also increases the maximum detectable Doppler frequency to 4 kHz, which can be further increased by phase-unwrapping. The input into the Kasai estimator can also be obtained from the Hilbert transformation as shown previously [5]. However, computing the Hilbert transformation requires 2n logðnÞ þ n complexity operations, where n is the number of samples per A-line, which is usually large ($4096) in order to digitize the carrier waveform. As a result, the computation time for a 200 Â 200 pixels image can be long ($7 s) [6], which is a potential limitation for real-time applications. In comparison, the time required for computing Eq. (2) is less than 0.08 sec for a 200 Â 200 pixels image. By shifting the I&Q data processing into analog hardware, significant time saving is achieved and real-time on-line signal processing is possible. Another approach is to use phase-locked-loop (PLL) techniques [24]. However, the minimum resolvable Doppler frequency is then large ($1 kHz), and the time required to achieve phaselocking is long ($80 ls), rendering this technique difficult to implement in high A-scan rate systems. In addition, these two parameters are an intrinsic trade-off pair in PLL. We conclude, therefore, that our current implementation is more sensitive, faster, and covers a larger velocity range for imaging microvasculature compared to existing techniques. Finally, the Kasai estimator could also be implemented in digital hardware [25] to further improve the computation speed.
Block tissue motion can be comparable or faster than the microcirculation velocities [26]. Fig. 4(a) shows an OCT structural cross-sectional image of human finger nail root region; the resulting de- gradation in the corresponding Doppler flow image is shown in Fig. 4(b). Note that the vertical bands indicate the block tissue movement velocity. In previous phase-resolved ODT systems [4][5][6], a glass window was pressed on the tissue to prevent such surface movement, which potentially alters the blood flow pattern. Previous retinal CDOCT [11] used a threshold technique to detect structural edges and a low-pass spatial filter to remove the motion artefacts. This method, like the cross-correlation registration method for motion artefact removal [16], depends on particular structural features and assumes such features are continuous across the entire OCT image, which may not generally be the case. Fig. 4(c) shows a typical velocity histogram, where the main peak indicates a block tissue motion velocity, and the much smaller sideband is due to blood flow of interest. Assuming the blood vessel diameter is much smaller than the A-scan depth, the number of pixels occupied by the moving blood would be less than that of the block tissue. In addition, the block tissue motion typically exhibits a narrow velocity distribution in the histogram, whereas the blood flow shows a wider distribution. The FWHM of the main peak in the histogram is $100 Hz, in agreement with the velocity resolution. To remove the block tissue motion artifact, the histogram is rotated within [Àp; p] such that the main peak is zero-centered, to produce a color flow image, as shown in Fig. 4(d). Our histogram-based technique exploits the excellent velocity sensitivity of phase-resolved ODT, assuming that small blood vessels are being imaged. It can remove motion artefacts both faster and slower than the blood flow velocity, provided that the block tissue velocity is fairly constant over the A-scans averaged (e.g., block tissue acceleration < velocity resolution/8-A-scan time, $3 cm=s 2 for our setup). Finally, the histogram method is not affected by aliasing due to the block tissue motion, and the estimated blood flow Doppler shift is unambiguous as long as the flow velocity does not introduce a phase-shift of more than AEp.
In summary, we have demonstrated the application of the Kasai velocity estimator in a phaseresolved ODT system, thereby improving the frame rate to approximately one frame per second for simultaneous structural and flow imaging in vivo. The minimum detectable Doppler frequency is $100 Hz, corresponding to a minimum detectable flow speed of $30 lm=s using eight axial-lines scanned at 8.05 kHz. The detectable Doppler frequency range is from )4 to 4 kHz, without using phase-unwrapping techniques. We have also implemented a histogram-based velocity segmentation technique to remove block tissue motion artefacts. Specifically, the improved ODT system allows measurement of not only the Doppler frequency but also the Doppler angle in vivo without physically immobilizing the tissue, and detected 1.4 mm/s blood flow in a 20 lm diameter capillary without motion artefacts.