One specific velocity color mapping using optical coherence tomography

Abstract. Depth resolved coherence gating along with Doppler shift detection of the carrier frequency is used for one predetermined velocity mapping in different flows. Bidirectional rapid scanning optical delay of optical coherence tomography system is applied in the reference arm. Tilted capillary entry is used as a hydrodynamic phantom to model a sign-variable flow with complex geometry. Structural and one specific velocity images are obtained from the scanning interferometer signal processing in the frequency domain using analog and digital filtering. A standard structural image is decomposed into three parts: stationary object, and positive and negative velocity distributions. The latter two show equivelocity maps of the flow. The final image is represented as the complexation of the three.


Introduction
Direction sensitive investigation of reversible flows using phase sensitive technique was applied in acoustic reflectometry 1 and optical coherence tomography (OCT). 2 . Velocity image reconstruction based on estimation of phase change between sequential A-scans is represented as the color or Doppler image in both. For that purpose, cross-3 and auto-correlation [4][5][6][7][8][9] functions of the reflected signal are used. For example, the OCT system presented in Ref. 3 is based on the calculation of cross-correlation function for a certain number, N, of adjacent A-scans. The mean phase shift variation is used as where «Arg» and «*»-denote the operation of taking argument of a complex signal and complex conjugation, respectively. S i;T and S Ã iþ1;T are the complex signals of adjacent A-scans. T is the temporal interval between the adjacent A-scans. The Doppler frequency is then represented as where f A ¼ ð1∕TÞ is the reference arm scanning frequency. Flow velocity, V, is then estimated as follows: where λ o is the central wavelength of the low coherence source, n is the refractive index of the media, and α is the Doppler angle.
The OCT system presented in Ref. 5 operates in the same way. The difference is that the cross-correlation function is used 5,6 in place of the autocorrelation Kasai function. 1 Flow velocity is then calculated from the phase shift of the consecutive A-scans. This allows obtaining high sensitivity to the motion of sharp surfaces or sparse particles in a flow, but leads to the 2π-ambiguity.
To reduce sensitivity to the changes of the incidence angle and to the velocity pulsations of the moving blood, Doppler spectra variance, σ 2 , can be used as 8,9 In Refs. 10 and 11, improved Doppler optical microangiography is described. In this method, the back reflected signal only from the moving particles is registered. This leads to a considerable reduction of the noise to compare with the standard backscattered signal processing. Blood flow velocity is determined from the phase shift of the signal carrier frequency.
A double correlation technique is notable for motion artifacts' reduction and background noise removal, 12 for which adaptive Wiener filtering is applied. It is also used to increase the structural image quality.
High sensitivity biological fluids' velocity measurement 13,14 is based on the phase synchronization of the reference signal, ΔΦ R , and in the sample, ΔΦ S . An electro-optic phase modulator with periodic phase tuning is used in the reference arm. 14 In practice, this time domain system allows the reduction of motion artifacts which causes interference pattern blurring, that, in turn, causes the sensitivity reduction of the flow registration. This approach is supposed to quantitatively register velocity change, but the Doppler frequency is, again, represented using the phase shift where k 0 is the central wave number and V R is the reference mirror velocity, which is defined as where ΔZ R is the reference mirror shift for the time, τ. To quantify ΔZ R , it is necessary to determine electrical values where β is the constant part of the phase, ΔV is the change of voltage for τ, and V π is the voltage at which the central wave number 14 is shifted by π.
Another approach to obtain velocity information is based on the speckle flickering frequency in the structural OCT images. 15 Speckle fluctuations depend on the average flow velocity. Therefore, the ratio of the high and low components of the intensity fluctuation power spectrum could show the average velocity of a flow. The Doppler shift is then calculated using the following equation: 16,17 where hi 2 AC i is the average value of the alternating photocurrent, hi 2 DC i is the average value of the direct photocurrent, and N S is the number of the detected speckles. And, it is assumed, the detector's surface should be uniformly illuminated.
The description method for one specific velocity (OSV) mapping is based on the registration of the Doppler shift of the carrier, f c , itself. This gives the average integrated information over the probe volume and does not depend on Doppler spectrum broadening. The latter does not give the absolute value of the flow velocity, but primarily has a statistical relationship with it rather than a functional one. In addition, OSV mapping does not have the 2π ambiguity disadvantage. The signal processing of Doppler spectra themselves is performed in the frequency domain after a short time Fourier transform (STFT).
Since it is difficult to register the Doppler shift of the carrier itself when applying ultrafast scanning techniques, 3-14 f D ¼ AE2V cosðαÞ∕λ 0 , we use the modified double pass rapid scanning optical delay (RSOD) 18,19 at lower scanning frequencies. RSOD allows the decoupling phase and group delay, as well as a shift of f c to the values of 20 to 30 kHz. The maximum signal-to-noise ratio of the detected photocurrent appears to be in the frequency range between the decreasing 1∕f noise and the increasing white noise. Using a time-domain optical delay and processing the signal in the frequency domain, it is possible to decompose a standard motionless structural image, and OSV images of a chosen predetermined velocity, V, with a certain predetermined accuracy, ΔV.
In the current work, the method of one predetermined velocity mapping is applied to study a flow with complex geometry in a tilted capillary entry (tilted die entry) 20 at a steeper angle, α ¼ 77 deg. The initial OCT image is decomposed into a structural image and two velocity images corresponding to the positive and negative directions of the flow. The final image is obtained after color coding and complexation of the three.

Materials and Methods
The experimental setup (Fig. 1) is a scanning fiber-based Michelson interferometer. Superluminescent diode (SLD) irradiation (λ o ¼ 1300 nm, Δλ ¼ 70 nm) after the fiber coupler goes to the reference arm-an optical delay line (modified RSOD 19 ). Another part of the irradiation goes to the sample arm for the investigated flow. The axial resolution of the system is about 11 μm. Tilted capillary entry is used to model a flow with complex geometry (Fig. 2). On the way back, irradiation is focused on the balanced detectors. After an amplifier and a tunable band-pass filter (2 to 200 kHz), the electric signal goes to a 10-bit analog-to-digital converter. Final processing is performed in the frequency domain using digital STFT.
For α ¼ 77 deg, the speed range in the hydrodynamic phantom is equal to 0 to 28 mm∕s. A 7.5-kHz Doppler shift corresponds to a velocity of V ¼ 17 mm∕ sec. The dynamic range of the system is about 100 dB. The system has been tested on the parabolic velocity profile in a cylindrical capillary with Doppler spectrum broadening Δf ∼ 0, 3-1, 6-kHz depending on the frequency of the centroid, f o The standard deviation of the parabolic profile ΔV∕V is about 5% to 7%. A hydrodynamic phantom has been made of a transparent plastic material with two different inner diameters d 1 ¼ 2.2 mm, d 2 ¼ 0.55 mm. This gives a complex flow with a converging ratio of 4∶1. The outer diameter of the phantom is 5 mm. This allows the study of the flow upstream (before the entry) and downstream (after the entry), with a resolution which cannot be detected by Doppler ultrasound systems. A 1% intralipid solution in distilled water was used for obtaining the backscattered signal and velocity images of the flow. The flow rate was equal to 0.11 ml∕ min, the average velocity of the fluid in the large lumen (d 1 ) was about 0.48 mm∕s and, assuming a parabolic velocity profile, the axial flow velocity was V 1 max ¼ 0.97 mm∕s. In the lumen with the smaller diameter (d 2 ), the average flow rate was about 7.7 ml∕s, and V 2 max ¼ 15.4 mm∕s. Since the percentage of the added intralipid was low, it was assumed that the density ρ and the dynamic viscosity of the solution were similar to those in water at 20°C.
The applied optical delay allows shifting the carrier to the region with minimal noise with a sufficiently fast repetition rate of A-scans (several kilohertz). This allows obtaining good enough structural images and registers the Doppler shift of the carrier corresponding to the indicated velocities. Odd A-scans always have a positive Doppler shift, while in the even A-scans, it is always negative. This approach gives the possibility of splitting the initial image into the three images corresponding to the steady object and images showing the positive and negative directions of the flow (positive and negative OSV images).

Velocity Mapping Algorithm
Using the aforementioned approach, one specific velocity mapping algorithm was suggested. It includes: (1) separation of the input data into two parts which correspond to the positive and negative shifts of the carrier frequency; (2) the independent construction of a two-dimensional (2-D) structural OCT image and (3) 2-D OSV images with (4) consecutive color coding.
x½s is the input set of data, where s ¼ 1;2; : : : ; S and S is the number of samples in the signal. It is received from the raw data file and represents the amplified and digitized mean intensity of the radiation detected by the system.
Then it is separated to the functions x odd ½p; k and x even ½p; k corresponding to the positive and the negative shifts of the carrier frequency as follows: x odd ½p; k ¼ x½p; k; for odd k 0; otherwise and x even ½p; k ¼ 0; for odd k reverseðx½p; kÞ; otherwise ; where p ¼ 1;2; : : : ; P, p is the count number within the A-scan, k ¼ S∕p, and P is the number of samples per one A-scan. The resulting signal x sum ½p; k is the algebraic sum of x odd ½p; k and x even ½p; k: x sum ½p; k ¼ x odd ½p; k þ x even ½p; k: Note that from the x sum ½p; k set the analog band-pass filter of the system separates the high-and low-parts of the spectrum.
To control the mapping process, two parameters have been chosen: the upper V upper and lower V lower limits of the mapped velocity V upper − V lower ≥ 0.5 mm∕c: On that basis, the upper ω upper and the lower ω lower cut-off frequencies of the band-pass filters were determined: They are necessary for determination of the frequency band from the functions x odd ½p; k and x even ½p; k.
To compare with the previously described algorithm, 20 this algorithm has better tuning of the chosen velocity diapason. It is determined not using the central velocity and its accuracy, but directly by independently choosing the upper and the lower velocity values. In addition, the ω upper and ω lower determinations have been also digitally improved. Currently, they are directly determined from the Doppler angle and the source central wavelength.
Functions x odd ½p; k, x even ½p; k, and x sum ½p; k are processed in three stages: (1) short-time Fourier transform; (2) Hilbert transform, and (3) taking the logarithm of the signal envelope. After that functions x odd ½m; k, x even ½m; k, and x sum ½m; k are obtained and considered. They represent three B-scans and correspond to the three parts of the image, Figs. 3(a) to 3(c) and 4(a) to 4(c), respectively. Note that m ¼ 1;2; : : : ; M, where M is the number of samples per one A-scan after the above processing.
Automatic identification of the directions of the flow x trends ½m; k is performed as the element-wise subtraction of x even ½m; k from x odd ½m; k with a shift by one column: x odd−1 ½m; k ¼ x odd ½m; k \ x odd ½m; k; x trends ½m; k ¼ x odd−1 ½m; k þ 1 − x even ½m; k; and subsequent logical analysis of the obtained results. Functions characterizing positively x pos ½m; k and negatively x neg ½m; k directed flows are determined as follows: x pos ½m; k ¼ 0; x trends ½m; k ≤ 0 x trends ½m; k; x trends ½m; k > 0 ; and x neg ½m; k ¼ x trends ½m; k; x trends ½m; k < 0 0; x trends ½m; k ≥ 0 : Then B-scans corresponding to the flows x pos ½m; k and x neg ½m; k are color coded considering the experimentally determined minimal level of the constant signal, L, which is separated from the noise (Figs. 3 and 4): x mpos ½m; k ¼ 0; x pos ½m; k ≤ L 1; x pos ½m; k > L ; and x mneg ½m; k ¼ −1; x neg ½m; k < −L 0; x neg ½m; k ≥ −L : After that functions x mpos ½m; k and x mneg ½m; k are complexed into a common array (OSV map) representing the flow direction and location x cart ½m; k [ Fig. 4(f)]: x cart ½m; k ¼ x mpos ½m; k þ x mneg ½m; k: Function x sum ½m; k is encoded to be at least 20 times higher than the functions x cart ½m; k and then is complexed with them: x last ½m; k ¼ x sum ½m; k þ x cart ½m; k: The resulting function x last ½m; k is normalized to the color scale, which includes shades of green to visualize the image of the phantom structure; red and blue colors are used to visualize the direction of the flow; the white color is used as a reference point (no movement, no reflection). To minimize complexation distortions, a proportionality factor between the x sum ½m; k and x cart ½m; k values equal to 20 was estimated as the best one. This allows simultaneously increasing the OSV accuracy values by 18% and the computation speed by 13%. The resulting image is presented in Figs. 3(h) and 4(g).

Results and Discussion
The described algorithm is realized in a form of a software package in graphical programming language «G» of LabVIEW. A simple model of a blood vessel was used for software testing and debugging. It is a transparent tube having an inner diameter of 0.55 mm. Figure 3 shows a structural image obtained by the system-(a); structural image of the tube-(b); structural image of a selected speed in the flow-(c) and final complexed image -(h).
The OSV image of the parabolic profile with other band-pass filters: 6 to 9 kHz (3.7 to 5.5 mm∕s)- Fig. 3(d

Conclusion
The method of Doppler mapping of a chosen velocity is applied on the basis of low coherence depth discrimination of the reflected signal and detection of the introduced carrier frequency changes. It gives 2-D maps with a quantitative functional (rather than statistical) relation between the Doppler shifts of the carrier and flow velocity without a 2π ambiguity. Additional color coding enables automated separation of the standard OCT structural image in to three parts: steady image and two velocity images, even in the more complex flow used to compare to the previously described image. 20 Note that velocity images are predetermined and quantitatively mapped and are proportional to the Doppler shift of the carrier itself. Colors are used as an additional option of the algorithm to separate a stationary object, as well as the positive and negative directions of the flow. It is possible to apply this approach in the study of small animals' hearts where similar flows with complex geometry appear in nature. Since color (Doppler) ultrasound systems have similar disadvantages as the low coherence optical ones, we suggest applying this algorithm in Doppler ultrasound diagnostics to study bigger biomedical objects.