Optical coherence tomography (OCT) was introduced more than two decades ago1 as a noninvasive modality for imaging transparent and translucent tissues with a resolution of a few micrometers.2,3 The first (and still dominating) application field of OCT was ophthalmology, where OCT revolutionized retinal imaging and diagnostics.4,5 In the later years, several functional extensions of OCT were developed, one of the most promising being Doppler OCT (DOCT)67.–8 which provides information on the movement of backscattering particles. Since major eye diseases such as diabetic retinopathy and glaucoma, as well as disorders such as retinal branch vein occlusion, are associated with alterations in blood perfusion, DOCT is of special interest for ophthalmic imaging.
A variety of different DOCT techniques have been reported in the literature. Examples are time domain-based optical Doppler tomography,9 phase-resolved DOCT (PR-DOCT),1011.12.–13 resonant Doppler flow imaging,14 joint spectral and time domain imaging,15,16 optical micro-angiography17,18 or single-pass volumetric bidirectional blood flow imaging.19
Apart from phase-based techniques, further principles can be used to obtain velocity information from OCT data. By using time varying speckle,20 it was possible to provide quantitative flow information for measurements on a flow phantom and in vivo. Although being able to detect motions normal to the optical axis, this technique cannot determine the direction of flow. Autocorrelation techniques21,22 utilize the statistical nature of intensity fluctuations of light backscattered by flowing particles. This technique was experimentally verified by using a flow phantom. Integrating dynamic light scattering and OCT (Ref. 23) is capable of distinguishing diffuse from translational motions and can measure axial and transverse components of the flow velocity. However, long measurement times are necessary and so far the computational effort for this technique is very high compared to conventional DOCT techniques.
In this work, we strive to obtain absolute quantitative and directional information on the flow, for which methods based on phase measurements are well suited. With these methods, the phase difference, , between signals of adjacent A-scans at partially overlapping positions is measured. The obtained is directly proportional to the axial velocity component, , at the measured sample position1) illustrates a limitation of common PR-DOCT: only the velocity component parallel to the direction of observation, , can be measured. This is indicated by the scalar product of the vectors and . If the angle between these two vectors is equal to 90 deg, the phase shift is zero and no velocity can be measured. Knowing can help to overcome the geometrical limitation and provides the absolute velocity
From the 3-D image geometry, the vessel orientation, and hence the angle , can be determined. However, if the measurements are performed in living tissue, motion artifacts during data acquisition frequently distort the 3-D geometry of the image, and the vessel orientation (and hence ) is unreliable.
Another option to overcome this problem is the use of two or more sampling beams simultaneously. Different variants of two- and multibeam setups have been reported for DOCT. Approaches that use two beams scanning the retina in parallel with a time delay were successfully demonstrated to extend the range of accessible velocities down to very low velocity values, thereby enabling improved contrast for small retinal vessels and capillaries.2627.–28 A combination of this approach with a statistical method based on temporal cross correlation was recently shown to provide quantitative and directional velocity information in flow phantoms.29
To solve the problem of the unknown Doppler angle directly, beams that sample the moving scatterers simultaneously from different directions can be used.3031.32.33.–34 Thereby, the individual components of the Doppler vector can be measured simultaneously and the angular dependence is avoided.
By using two orthogonally polarized beams produced by a Wollaston prism, the first dual beam DOCT setup with bidirectional probing was realized.30 Successful measurements were demonstrated for the flow profile in glass capillaries. Using a trigonometric approximation in order to account for the orientation between illuminating beams and capillaries, the system was independent of the angle of incidence for angles close to 90 deg. One drawback of this approach was that only two beams were used, providing just two components of the velocity vector, therefore, additional information on the capillary orientation was necessary to obtain the absolute velocity. This technique was further improved and in-vivo applications in the human eye were reported.34,35 Although still limited to measure only two components of the velocity vector directly, the third dimension could be retrieved by measuring the azimuth angle of the retinal vessel from a fundus image. Recently, this technique was extended by rotating the orientation of the plane spanned by the dual-beam geometry with a mechanically driven Dove prism.36 Proper alignment of the beam geometry with respect to the imaged blood vessels can provide the absolute velocity vector from just two sampling beams. However, the proper alignment of the Dove prism complicates the measurement since it requires a priori knowledge of the investigated vessel geometry.
PR-DOCT systems acquiring all the three velocity vector components, without using additional structural data for the velocity reconstruction, have been reported on a very limited basis in the literature.31,33
A very compact design used a single beam and divided its cross section into three segments with a beam divider consisting of a glass plate with sectors of different thickness. After directing the different cross-sectional segments of the beam through a focusing lens, different measurement orientations were achieved. The velocity information corresponding to the three components is encoded in depth due to the different optical path lengths of the components within the three-segment beam divider. This approach showed promising results for in-vitro applications.33 An advantage of this method is that the three velocity components along one depth profile can be recorded with a single A-scan. This, however, severely limits the achievable depth range.
The purpose of this work is to present a new improved PR-DOCT system that overcomes the limitations mentioned above. The instrument uses three separate sampling beams which provide access to the three components of the velocity vector. This enables the reconstruction of the absolute velocity vector without requiring any additional information on vessel orientation from structural data, thus avoiding any artifacts due to the sample motion. Furthermore, the use of three independent channels in parallel effectively eliminates any cross talk37 and provides full ranging depth for each channel. We demonstrate the performance of the system by calibration measurements in a rotating disk, by measurements in a capillary flow phantom of well-defined flow, and by analyzing the flow in a venous bifurcation in the retina of a healthy volunteer.
Sampling beam geometry
Our three-beam DOCT system employs a special sampling beam geometry depicted in Fig. 1. Three parallel sampling beams, separated from each other by equal distances, illuminate a joint focusing lens which focuses them to a mutual spot on the sample. This yields the geometry of an equilateral triangle-based pyramid, with the tip of the pyramid being the mutual focal spot of the respective beams (see Fig. 1). Knowledge of the distance between the parallel beams and the focal length of the lens provides all information necessary to determine the exact spatial orientation of the three measurement beams. Reconsidering Eq. (1) and expressing it for each measuring direction yields a system of three equations with three unknown coefficients:
The unknown coefficients in Eq. (3) are the three components of the velocity vector, , , , and can be determined by solving the system of equations.
Our system is based on spectral domain (SD)-OCT (Refs. 3839.–40). A schematic drawing of the three-beam DOCT system is shown in Fig. 2 and a description can be found in the following. The experimental setup combines three SD-OCT systems. We employ three similar superluminescent diodes (SLD) (EXS0840, Exalos, Switzerland) emitting light centered at a wavelength around 840 nm with a bandwidth of 50 nm. The SLDs are coupled to single mode fibers, guiding the light to a set of miniature collimators (beam diameter ). An in-house designed mount, machined out of one piece of brass, aligns the laser beams according to the corners of an equilateral triangle. This mono-block construction avoids temperature and vibration drifts between the beams. The beams are aligned parallel to each other with a beam distance of 4.4 mm. The three beams are directed toward a nonpolarizing bulk beam splitter which separates the respective beams into the sample and reference arms. For in-vitro measurements, we use a (R/T) splitter, while for in-vivo applications in the eye, we use a (R/T) splitting ratio.
In the sample arm, two telescopes are used in order to set the beam distance. The first telescope, , has a focal length ratio of and is used to achieve a narrow beam spacing at the galvo scanner. Having three beams located at the corners of a triangle, it is impossible to provide all beams to be at the pivot axis of the scanner. An off-axis position introduces phase shifts due to the scanner movement.41 Decreasing the beam distance simultaneously reduces the phase shift and minimizes fringe washout. Remaining phase shifts can then be compensated for in postprocessing.
With the second telescope, , featuring a focal length ratio of , the beam distance is increased to a value satisfying the needs for in-vitro and in-vivo measurements. A distance of 3.5 mm provides a substantially different inclination of the measurement beams with respect to each other and simultaneously allows for the penetration of a dilated human pupil with all three measurement beams.
For the in-vitro case, an objective lens focuses the beams to a mutual spot on the sample while for in-vivo measurements this is performed by the refractive optics of the human eye. It is important that the optical axis of the final lens and the equilateral triangular-based pyramid height coincide to produce the measurement geometry described above. The distance between the parallel beams after the second telescope, , and the focal length of the final lens, or the eye length, determines the actual orientation of each measurement direction. In order to ensure the correct measurement geometry for in-vitro measurements, the beam distance was checked by a charge coupled device (CCD) camera at several positions before and after the mutual focusing lens. The obtained actual geometry was in very good agreement with the expected values.
In the reference arm all three beams pass a neutral density filter, which is used to set the light power at the detection unit close to the saturation level of the CCD line-scan cameras.
Finally, backscattered light from the sample is recombined with the light reflected at the reference mirror at the beam splitter and coupled to a second set of miniature collimators. Single mode fibers guide the light toward three identical spectrometer units, consisting of a collimator, blazed reflective grating (), focusing lens (), and a line-scan camera (Atmel Aviiva, 2048 pixels). During all our measurements, the cameras are operated at a line-scan frequency of 20 kHz.
The system combines three individual SD-OCT setups and requires an accordingly higher effort in terms of setup building and alignment. Especially, a careful alignment of the spectrometer units to provide similar phase characteristics with depth of the Fourier transformed signals in all three channels is necessary. The final bulk size of the entire setup fitted on an average sized optical table () as can be seen in Fig. 3. However, while the system is presently rather bulky, the complexity may be reduced in future versions of the instrument by replacing the three separate detection units by a single camera system with joint spectrometer optics.4243.–44
For the in-vitro measurements an illumination power of 700 μW per beam is set, adding to a total illumination power of 2.1 mW and leading to a measured sensitivity (close to zero delay) of 98 dB per channel. However, for in-vivo measurements, the illumination power is decreased to 700 μW in total or 233 μW per beam, which is below the lasers safety limits.45,46 With the lower illumination power, a sensitivity of 95 dB was measured for each channel.
All three beams measure the same spot at the sample simultaneously. During postprocessing, the three measurement directions can be regarded as independent and each channel is first evaluated individually in terms of intensity and Doppler data. Then, the velocity data along the different measurement directions is combined to compute the absolute velocity vector according to Eq. (3).
In a first step, intensity images are generated by standard SD-OCT preprocessing procedures (mean spectrum subtraction, rescaling from wavelength to wavenumber space, zero padding, dispersion correction, and inverse Fourier transform). In a second step, the phase difference of adjacent A-scans is calculated. In order to eliminate background phase noise, we compute this phase difference only for pixels above a certain intensity threshold level (typically four times the intensity noise level). This leaves us with three independent phase difference tomograms of the same position in the sample. Due to the mono-block construction of the miniature collimator mount and the bulk optics interferometer setup, our system is very phase stable. No phase drift correction is needed along one B-scan. However, it is necessary to correct for the phase shift introduced by the galvo scanner. To do so, we apply a histogram-based method13 for each channel. Under the assumption that more static than dynamic tissue is imaged, we compute a phase histogram of an entire B-scan. This histogram reveals the phase shift introduced by the scanner and bulk motion (e.g., head movement during the measurement). Subtracting the most populated histogram value from all the other pixel values in the B-scan results in a phase shift corrected image.
To improve the signal quality, several B-scans are averaged (for each channel separately) by an algorithm working in complex space.47 Before the averaging, B-scans have to be aligned to eliminate bulk motion artifacts. For in-vitro measurements, no sample bulk motion is present, making a motion correction unnecessary. The averaging is typically applied over 20 B-scans. For in-vivo evaluation, a basic eye motion correction is performed, which correlates the position of vessel centers within the different B-scans to each other. Using the number of B-scans that are recorded during the time of one heartbeat, we average over the systolic and diastolic phase of one pulse. This typically led to 15 averaged B-scans for the in-vivo measurement.
The index runs from 1 to 3 and indicates the individual measurement directions while the indices and indicate the position of the respective pixel in the B-scan. According to Eq. (4), the averaged phase images are rescaled to the velocity measured along the respective beam direction for every pixel in one B-scan. This yields the axially measured velocity components , , .
By using the known beam geometry (vectors , , ), we can now solve the system of equations [Eq. (5)] for each pixel position with its individual set of axial velocities, . This leads to the distribution of the three velocity components, , , of the flowing particles. The total flow can then be obtained by integrating the velocity values over the scanned area of interest. The postprocessing algorithm is outlined in a flow chart in Fig. 4.
For postprocessing, a PC featuring an Intel i7 CPU (3.2 GHz) with 24 GB RAM was used running a 64 bit Windows 7 operating system. Postprocessing of a typical data set, consisting of 50 B-scans (1024 A-scans each), with a nonoptimized Labview program took approximately 5 min.
To demonstrate the performance of our three-beam DOCT system in vitro, we performed measurements on a rotating disc and a flow phantom. In a next step, we conducted in-vivo measurements on human retinal blood vessels.
In order to prove the functionality of our setup and to calibrate the instrument in a simple and well defined experiment, measurements on a rotating and tiltable disc were performed (cf. Fig. 5). The stationary movement of the disc made it possible to avoid scanning.
For the computation of the velocity vector components, knowledge of the exact beam geometry plays a crucial role in our approach. We used the following calibration procedure to ensure correct beam geometry: The galvo scanners were set at a fixed position and a constant angular velocity was set for the disc. By choosing a measurement spot at a certain radius, the velocity vector was defined. Knowledge of the velocity vector provides the expected phase differences for the respective measurement directions (channels) during that measurement. In a next step, the focusing lens was aligned orthogonal to the optical axis while the values were monitored. After the expected phase differences for all three measurement directions were achieved, the optical axis and the equilateral triangular pyramid height were assumed to be collinear and the measurement geometry was regarded as correct.
After the alignment procedure, measurements were performed for various settings. At each setting, 1024 A-scans were recorded and phase differences were averaged in the complex space. The averaged values of were used to calculate the three axial velocity components corresponding to the three measurement directions [Eq. (1)], and with equation system (3) the velocity components , , were obtained.
For the first test, the disc was aligned orthogonal to the optical axis (i.e., parallel to the plane) and the angular velocity was kept constant. At a distance of 9 mm from the disc center, points separated by 45 deg were measured. Figure 5 shows a sketch of the alignment of beams and rotating disc as well as an overview regarding the measurement points.
Figure 6 shows the results. Velocity components , , and the absolute velocity are shown as a function of sampling point position along a circle with the . The preset absolute velocity for this angular velocity and radius was . This is in good agreement with the measured absolute velocity of (). As expected, and oscillate in a sinusoidal manner along the circle circumference and are 90 deg out-of-phase while is close to zero. For estimating the error of , , and velocity components (error bars in Fig. 6), we took the standard deviation for the absolute velocity () as the highest expected error for the respective components.
In order to demonstrate the performance of our method to measure the velocity component, we conducted a further test where the disc was tilted between and around the -axis in steps of 2 deg. Measurements were taken at a distance of 3 mm to the center of the disc. A schematic of the test setup and the obtained data can be seen in Fig. 7. The preset absolute velocity at the measurement position was . The measured absolute velocity of is in good agreement. Velocity components , , and as well as the absolute velocity are shown as a function of tilt angle . As expected, shows an approximately linear increase with within this range of the tilt angle, while slightly decreases and remains close to zero (the slight increase of with is probably caused by slight deviations of the measured position from the -axis). For an approximate error estimation, we used again the standard deviation of the absolute velocity.
To demonstrate the capability of our method to measure flow quantitatively, we performed measurements in a flow phantom consisting of a glass capillary perfused with a scattering fluid. To illustrate the velocity vector reconstruction for all three-dimensions (3-D), the capillary was rotated in planes orthogonal and parallel to the optical axis.
Outer and inner diameters of the capillary were 1.0 and 0.3 mm, respectively. As a scattering medium, we used milk diluted with water in a ratio of . A special mount, providing three rotational degrees of freedom, was used to align the sample. In order to maintain an accurate and constant flow in the capillary during the measurement, an injection pump for medical applications was used (MGVG Combimat, adjustable flow range leading to mean velocities of for a capillary diameter of 0.3 mm).
Operating the system at a line scan rate of 20 kHz, we recorded 20 B-scans consisting of 1024 A-scans at a constant B-scan position. The scanned length was 1.2 mm. Scanning changes the beam geometry as the measurement beams move through the optics and the sample. These geometry alterations can lead to systematic errors in the velocity calculation. We considered these effects by calculating the beam geometry as a function of galvo scanner position. This procedure yields new beam geometry vectors (, , ) for each measurement spot on the sample. These vectors were then inserted into our standard postprocessing procedure. After compensating the effect in postprocessing, we compared the obtained results with the original values. No significant difference was observed for the small scan angle used (). Therefore, we omitted the described beam geometry correction for small scanning angles.
For the measurements, the mutual focal plane of the three sampling beams was aligned in the center of the capillary while the zero-delay was shifted to the capillary edge. The calibration procedure was similar to the one performed for the rotating disc. For a known capillary orientation, a constant flow was set at the injection pump. Without scanning, the signal phase shift obtained from the scattering fluid along a single depth profile was monitored. The capillary orientation defined the velocity vector orientation and the focusing lens was aligned until the expected phase differences for all three channels were achieved.
In a first experiment, we demonstrate the velocity vector reconstruction within a plane orthogonal to the optical axis. Figure 8 shows a schematic of sample and beam geometry. A constant flow velocity of (mean velocity over the whole profile) was set at the injection pump. Measurements were taken at different settings of the capillary orientation angle , varying between 20 deg and 160 deg. A set of images illustrating the reconstruction method for a measurement at can be seen in Fig. 9. The left-hand column [Fig. 9(a)] shows intensity B-scans obtained in the three channels. The next column [Fig. 9(b)] shows the corresponding averaged phase difference images of the capillary center. Gray pixels indicate areas where the intensity threshold was not met (essentially indicating the glass walls of the capillary). The third column [Fig. 9(c)] shows the reconstructed velocity components , , . Because the capillary was oriented at , the components and are equal, while (since the capillary is in the -plane), as expected. Figure 9(d) shows a 3-D plot of the absolute velocity , showing the expected parabolic velocity profile.
Figure 10 shows a plot of (mean) velocity components (blue), (green), (orange) and of the absolute velocity (red) as a function of angle . The expected values are indicated by lines, the measured values by circle symbols. The error bars equal the standard deviation of the absolute velocity, calculated over the data obtained at the different settings of . As can be seen, there is some discrepancy between expected and measured values (on average 15% for the absolute velocity). One reason for this discrepancy might be the refraction of the sampling beams at the glass capillary which causes a local deviation from the expected beam geometry.
In a second experiment, we demonstrate the velocity vector reconstruction within a plane parallel to the optical axis (-plane). The orientation angle was varied between 0 deg and 30 deg. The flow within the capillary was set to . A sketch of the beam and sample geometry can be seen in Fig. 11.
Figure 12 shows the results for a capillary orientation of . All three phase difference tomograms, , show a large phase shift [Fig. 12(b)]. This can be explained by the fact that the velocity components along the directions of the three sampling beams are now considerably larger than for the case of the capillary within the -plane. As expected, the reconstructed velocity component is close to zero, while and have nonzero values [Fig. 12(c)]. A nearly parabolic velocity profile of the absolute velocity is again visible [Fig. 12(d)].
Figure 13 shows a plot of velocity component (blue), (green), (orange) and of the absolute velocity (red) as a function of tilt angle . Lines show expected values, circle symbols measured values; the error bars equal the standard deviation obtained for the measured absolute velocity values across all settings of . In comparison to the preset absolute velocity of , the measured value of shows a good agreement. Deviations between expected and measured values are again likely caused by a distortion of the beam geometry by refraction at the glass capillary.
To demonstrate the functionality of our setup for in vivo retinal imaging, measurements were performed in the eye of a healthy volunteer after full informed consent was obtained. The axial eye length was measured by partial coherence interferometry48 (Zeiss IOL Master) to calculate the exact beam geometry for the flow measurement. To ensure the correct beam geometry, the focal plane position was adjusted by shifting the second lens of the second telescope (cf. Fig. 2) along the optical axis. A correct setting of the mutual focal plane of the three beams was achieved when the volunteer confirmed the fusion of three scanning lines to a single line. Furthermore, an online preview was used in order to ensure that the respective blood vessels were in the same position for all three measurement channels.
Imaging was done at an A-line frequency of 20 kHz. The scanning distance was approximately 1 mm or 3 deg. For one measurement, we recorded 50 B-scans at the same position comprising 1024 A-scans each. The subject heart beat rate was determined to be during the measurement, and in the postprocessing an artifact free set of 15 consecutive B-scans, covering the time of one heart beat cycle, was chosen for further evaluation. Before averaging, B-scans were aligned to correct for eye motions.
Figure 14 shows a fundus image of the measured eye. The aim was to quantify the total in- and outflow at the indicated venous bifurcation. For that purpose, measurements were performed at positions , , and . A typical set of tomograms obtained by the measurement at the venous inflow position can be seen in Fig. 15. The second column [Fig. 15(b)] shows averaged phase difference images of 15 B-scans. The vessel can be easily distinguished from the static tissue in all three channels. The velocity component images [Fig. 15(c)] show larger velocities only for components and , while is close to zero, as expected for a vessel away from the optic nerve head. and were averaged over the vessel cross section and the velocity vector in the -plane was determined. The vector angle of 135 deg to the -direction at is in good agreement with the vessel orientation at this point in the fundus photo. Flow rates , , at , , and were calculated by integrating the velocity components over the vessel areas.
Figure 16 illustrates flow orientation and quantity for the three vessel branches of the bifurcation obtained at the three measurement locations. The length of the respective arrow is proportional to the flow quantity (in ); the arrow orientation corresponds to the retrieved flow vector orientation. The vessel diameters were measured to be between 0.12 and 0.17 mm while the total flow rate in the respective vessels showed a range of .
These results are in good agreement with the results found in the literature: measurements done with laser Doppler velocimetry (LDV)49 suggest total flow rates of for the same vessel diameter range. A more recently published paper24 using a Fourier domain D-OCT system reported flow rates of for vessel diameters of 0.070–0.152 mm.
Table 1 summarizes the results. The total measured inflow () is in good agreement with the total measured outflow ().
Summary of blood flow measurement results on a human retinal venous bifurcation.
|Inflow Q1||Inflow Q2||Outflow Q3|
Discussion and Conclusion
We developed a new three-beam DOCT system that allows the 3-D reconstruction of flow orientation and the determination of the absolute flow velocity, without the need for information on the orientation of the vessel. We demonstrated the functionality of the new system in a simple test object, a glass capillary, and—to the best of our knowledge for the first time—in the retina of a human volunteer in vivo.
By measurements performed on a rotating disc, we demonstrated the functionality of our system for velocity vector reconstruction in all 3-D. Scanning a single point on an evenly moving object, a good correlation between preset velocities and measured velocities was observed. It should be mentioned that with this simple test object, no scanning was performed. In the more realistic flow phantom, we tested the performance of our system in a real imaging situation, i.e., with beam scanning in operation. By pumping a scattering fluid through the capillary, a constant flow was provided for measurements. Aligning the flow in a plane orthogonal and parallel to the optical axis, we showed that our approach allows one to reconstruct the velocity vector in all 3-D for arbitrary velocity orientations. The observed deviations between expected and measured flow values are probably caused by refractions at the glass capillary which distort the sampling beam geometry. Since the refractive index differences between vessel walls and surrounding tissue are considerably smaller than those between the capillary glass and its adjacent media (air and aqueous solution) we expect that this effect will be lower for measurements in tissue. However, since the local distribution of refractive index in tissue is not exactly known, this effect might cause minor deviations of velocity measurements from true values.
Applying the new three-beam DOCT technique to the human eye yielded retinal blood velocity and blood flow results which are in the same range as those measured by other methods.24,35,49,50 Furthermore, the good agreement between bifurcation in- and outflow is very promising and supports the reliability of the measurements. Present shortcomings comprise the rather low transversal resolution of , caused by the low numerical aperture of the beams used, and the presently complicated alignment procedure for in-vivo measurements that requires feedback from the subject. The former problem might be mitigated by using collimators with a larger beam diameter (although there is a trade-off with beam separation), and the latter problem might be solved by adding a fundus view capable of resolving the laser beams on the retina for quick and objective focusing.
To summarize, we developed a new three-beam SD-DOCT system and demonstrated its ability to reconstruct the velocity vector in 3-D. The instrument was demonstrated in test samples and in the retina of a human volunteer in vivo. Reliable quantification of the flow was demonstrated by measuring in- and outflow at a venous bifurcation.
Financial support from the European Union (project FUN OCT, FP7 HEALTH, Contract No. 201880) and from the Austrian Science Fund (FWF Grant No. P19624-B02) is gratefully acknowledged.