|Home | About | Journals | Submit | Contact Us | Français|
Doppler optical coherence tomography has the capability to measure blood flow quantitatively and in vivo. As only the axial component of the velocity can be assessed, the measurements have to be corrected for the angle of the vessels. We present a novel approach to extract quantitative flow data from circumpapillary scans in vivo on the human retina by registering the circular scan to a reference volume scan and extracting the angle directly from the volume. In addition, we perform phase unwrapping and interpolation of the flow under the assumption of a parabolic flow profile. We demonstrate the repeatability of the methods by applying it to different retinal vessels, achieving coefficients of variation of the average velocity of 3 to 8%. Results on the pulsatility and resistance index are also presented.
Optical coherence tomography (OCT) has already become an established diagnosis tool in ophthalmology . Standard OCT imaging uses backscattered intensity as contrast mechanism for displaying retinal structure. Fourier domain OCT brought the capability to image full 3D sections with high resolution and high morphological reliability. Outstanding sensitivity allows for imaging speeds to sample full volumes in fractions of a second [2,3]. However, intensity contrast has its limitations for displaying tissue physiology, which is important for understanding the pathogenesis of important retinal diseases as well as for following treatment progression. Recently, functional extensions of OCT such as polarization contrast OCT, and Doppler OCT (DOCT) came into focus of research . DOCT has the potential to diagnose retinal diseases at an early stage, since perfusion is highly sensitive to any disorders. The latter has been demonstrated with Laser Dopper velocimetry for various diseases [5–7]. DOCT does not only serve as contrast tool for tissue vascularization down to the capillary level [8–12], but yields also full quantitative information on retinal perfusion within 3D volumes, i.e. quantitative 4D angiography [13–15]. Sharing the high-speed advantage of FDOCT [16,17] it is possible to follow flow dynamics over time within individual retinal vessels [18–20].
The measured Doppler shift is proportional to the flow velocity and their relation is given by the formula 
where is the group refractive index of the medium, v is the absolute flow velocity, α is the angle between the OCT beam and the flow, and is the central wavelength of the light source. The dependence on the angle α makes it difficult to extract the correct flow information. As pointed out earlier, the correction of blood flow information and the extraction of quantitative blood flow parameters is very important as it may have a huge impact on retinal diagnosis. Presented solutions of this problem include the direct reconstruction of vessel orientation from 3D volumes [22,23] or the usage of multi-beam illumination to obtain the 3D velocity vector [24–28]. Less complex methods record two arc scans with increasing radius or use a double-circular scan pattern in order to extract the local vessel gradient from adjacent vessel cross sections [29,30]. Those methods rely on the exact knowledge of the distance between the double scans for extracting the angle α. However, during in-vivo imaging, any patient motion will distort the scanning pattern, and thus impair the correct angle determination. Still, the advantage is the possibility to assess perfusion dynamics within vessels, since circular or arc scans allow for good temporal sampling of the pulse cycle. Also, circumpapillary scans yield important information on the total retinal arterial in- and venous outflow.
For tackling motion artifacts, we therefore propose to register circular scans to a reference OCT volume recorded just prior to the circular scans. The vessel gradients can then directly be read off the 3D volume by establishing a look-up table for angles of individual vessels within the volume. The 2D-3D registration should be virtually free of motion artifacts, since the exact position of the circular scan can be retrieved.
The optical setup is shown in Fig. 1 . The light of a super luminescence diode (SLD) with 65 nm spectral width, centered at 860 nm (Exalos Inc.), resulting in 5 µm axial resolution (in air), is fiber-coupled into a bulk-optics Michelson interferometer. In the sample arm, the light is directed to the eye after passing a X-Y-galvo-scanner unit (SC). The diameter of the collimated beam at the cornea is 1.8 mm resulting in a spot size on the retina of 14.7 µm. The power of light at the cornea was measured to be 610 µW, which is at this wavelength well below laser safety standards. The reference arm contains a dispersion compensating prism-pair (DISP), a water chamber (25mm), and a mirror mounted on a linear translation stage (M) which is used only once to adjust the coherence gate for retinal imaging.. At the exit of the interferometer the light of the two arms is recombined and guided to the spectrometer through a single mode fiber. The spectrometer consists of a volume phase holographic transmission grating (Wasatch, 1200 lines/mm, denoted by DG in Fig. 1), a camera lens (f = 85 mm, denoted by L in Fig. 1.), and a high speed CMOS line scan camera (Basler Sprint spL4096-140km, denoted by LSC in Fig. 1.), with 4096 pixels. We imaged the spectrum of the SLD onto 1152 pixels, allowing for an effective maximum camera speed of 200 kHz. All measurements presented in this paper were however acquired at a line rate of 30 kHz. This line rate has been chosen in order to match the dynamic range of detectable flow speeds with the flow speeds in the particular retinal region . The recorded volumes cover a patch on the retina of 6.45mm x 6.45mm x 2.3mm (1000x100x575 pixels). The recording time of a volume was roughly 3 seconds. The circumpapillary scans have a diameter of 4.2 mm and contain 3000 A-scans. This results in a lateral oversampling factors of 2.27 for the volume scans and 3 for the circular scans.
All experiments were performed in compliance with the relevant laws and institutional guidelines and were approved by the ethics commission of the Medical University of Vienna.
Before registering the circular scans to the volume, we determine the vessel orientations within the volume. Like this, we can create a look up table for the angles of the vessels at their respective locations. Clearly, the proper registration of the recorded B-scans is important for generating a reference volume exhibiting physiologically correct vessel orientations. For this task we use two fast reference scans in Y-direction  prior to raster scanning as shown in Fig. 2 . The fast reference tomograms can be regarded as practically free of motion artifacts. They consist of 1000 A-scans and are recorded at 30 kHz within only 0.033 seconds. Their size is 20% larger than that of the raster-scanned volume in this direction for registration purposes. We then need to find the position of the subsequently recorded volume with respect to the reference tomograms. For this purpose, we extract virtual tomograms in Y-direction from the volume and correlate the intensity-projection along the Z-axis with that of the respective reference tomogram. In order to quantify the agreement between both tomograms, we determine the correlation coefficient between respective A-scans and calculate its mean value across the tomograms. The scan with the maximum mean correlation is chosen as the tomogram, which correctly represents the reference scan. This is repeated for the second reference tomogram. Since the fast scanning direction of the raster scan can be regarded as free of motion artifacts we can align the height of the left and right side of each volume B-scan along the reference tomograms in Y-direction. The final registered volume then has the correct orientation with respect to both reference tomograms (see Fig. 2).
So far only the intensity information has been used for registration. The segmentation of vessels is done with the help of the actual Doppler tomograms which are obtained by calculation of phase differences between adjacent A-scans using a circular averaging technique . Before segmenting the vessels, we need to correct the phase-difference image for bulk motion artifacts using a histogram-based procedure [3,9]. We then manually segment selected vessels from the volume scan starting on a small circumpapillary circle and moving outwards at minimum step size [see Fig. 3(a) ]. In order to achieve this an ellipse is drawn to describe the vessel cross-section, the center of which is chosen to be the center of the vessel.
The connected center positions describe the vessel traces through the volume. After smoothing the traces by applying a running average to each trace-projection, we can determine the local vessel orientation α with the following formula:
The distances Δx, Δy and Δz are defined by the volume sampling in the respective dimensions. The angles are finally stored together with the corresponding positions in a look-up table for each vessel trace.
The different sampling in X (B-scan)- and Y-direction has to be taken into account when extracting virtual scans from the volume, be it in this step for the angle extraction or at a later stage to register the circular scans to the volume. The difference in sampling results in virtual circular scans having worse resolution on the sides that are in Y-direction (nasal and temporal) as compared to the sides in X-direction (inferior and superior). Hence, nasally and temporally, the A-scans are repeated for different angles as a consequence of the reduced sampling. However, the characteristic vessel shadows are little affected by this unequal sampling which is also why it did not have a noticeable effect on the registration.
The different sampling of the volume along the X- and Y-axis was also taken into account by transforming the differences in [pixels] to real distances in [mm] before applying Eq. (2).
In order to quantify the uncertainty of our manual segmentation, we have repeated the manual segmentation of two whole vessels four times - two experts independently segmented two vessels two times each on 10 different virtual scans with increasing radius. The mean standard deviation of the angle segmentations over the two vessels was found to be 1 degree.
The central idea of this work is to register the recorded circular scans to a reference volume and extract the previously calculated angles at the respective positions [see Fig. 3(b) and 4(a)
]. For this task we use intensity projections, since the enhanced scattering in blood leads to shadowing effects visible as dips in the projected intensity trace across a tomogram [cf.
Fig. 4(b)]. Those distinct features within such a shadowgram can be used to guide the registration. As a similarity measure, we use the correlation coefficient between the shadowgrams of the recorded circular and extracted virtual scan. The registration is facilitated by not projecting the whole B-scan structure along the Z-axis but only the photoreceptor layer (PR) and retinal pigment epithelium (RPE) [indicated by red lines in Fig. 4(c)]. These layers are segmented by a simple thresholding procedure. For stable registration results the shadowgram should exhibit strong variations in the area of the vessels and small variations in areas of bulk tissue. The vessel shadow features can be enhanced by Fourier band-pass filtering of the shadowgram. So, our cost function becomes
where f and g denote the projection of the circular and the virtual scan, respectively, and h stands for the band pass filter applied to the signals, i.e. where the asterisk '' stands for discrete convolution.
The correlation coefficient r(.) is defined as
where denote averages.
Clearly, we cannot a priori assume that the circular scanning pattern will be imaged without deformation onto the retina. With good approximation we can consider an arbitrary elliptical pattern. As parameters in the optimization we use the center point of the ellipse and the radii in X- and Y-direction. We found that the rotation angle of the ellipse about its center did not affect the outcome of the registration in the present experiments. The dependence of the cost function on the parameters in the optimization is hidden in the functions f and g. We use the Nelder-Meads simplex algorithm to minimize the cost function , which in our experiments gave much better results than Newton-based methods. The motivation is that we have a multi-dimensional optimization landscape with a nonlinear dependence of the cost function on the parameters that are varied during the optimization. The Nelder-Meads optimization algorithm is based on the concept of a simplex the edges of which are reflected, contracted, or expanded in accordance to certain heuristic rules. This makes the simplex wander along the optimization landscape towards the optimum. It shall be noted that the simplex used in the Nelder-Meads algorithm has nothing to do with the well-known Simplex algorithm from linear programming. Figure 3(b) shows a circular Doppler scan on top of the fundus image in order to illustrate the extraction of the angles after the registration. The red lines indicate the vessel centers that have been determined. An example result of the procedure can be seen in Fig. 4. In our experiments this yielded a very stable optimization for finding the best fit between actual circular and extracted virtual scan. The optimization time on an Intel Core i7 CPU 920 system with 2.67 GHz was between 1 and 2 minutes per slice.
Once the correct position of the actual scan within the reference volume is found, one can directly assign the angles to the vessel cross sections of the scan. The registration runs fully automated for any number of circular scans and no intervention is necessary. After having registered circular scans to the volume, we can select specific vessels for studying their dynamic flow characteristics over time. This is done via choosing the respective minimum in the projection of the PR and RPE [see Fig. 4(b)] for the vessels that we want to evaluate.
In order to ease the quantitative extraction of perfusion parameters from vessel cross sections we need to register the vessels along time.
As a result of motion artifacts the lateral distance between vessel cross sections as well as their position within circular scans might change between successive scans. This is best seen from the intensity projection of the circular scan series over time [Fig. 5(a) ]. The dark bands in Fig. 5(a) correspond to vessel shadows. This representation allows for manually choosing an area of interest for which the respective vessel cross-section will be registered over time using cross-correlation. An example of the registration output can be found in Fig. 5(b).
After the registration, the elliptic area of interest of the vessel needs to be chosen only once and the flow analysis can be performed automatically over time. Figure 6 (Media 1) shows a screenshot of a video that illustrates the different steps, also showing the segmentation of the flow by using ellipses.
Large flow values especially in arteries cause the appearance of phase wrapping artifacts when the flow value exceeds the unambiguous velocity range , where T is the time between successive A-scans and ng denotes the group refractive index of blood which we have assumed to be 1.37 (the value was calculated from reference , p. 141, as an average of the values given for 633 nm and 1080 nm).The calculated phase difference ΔΦ between adjacent A-scans is mapped to the range [-π, π]. Single phase-wraps can be easily corrected for by determining the axial flow direction from the Doppler value close to the border of the vessel. If we know that the actual phase difference , we have
in case of , we have
After phase unwrapping the phase shifts are mapped to the range [0, 2π] or [-2π, 0], respectively, depending on the direction of the flow. Double phase wraps did not occur during our experiments. Another artifact for Doppler OCT using spectrometer based systems is fringe washout. This effect means the movement of interference fringes across the detector during exposure and thus OCT signal loss for moving structures. Fringe washout manifests itself as missing information especially at vessel centers where flow values are large [Fig. 7(b) ]. This artifact is corrected for by the assumption of a parabolic flow profile and least squares interpolation of missing Doppler data. This leads us to the following set of equations:
where a stands for the parameters of the paraboloid that need to be fitted and the coordinates of the measured points are denoted by (x,y,z). This system of equations can be solved in a least-squares sense.
The main steps of the pipeline are presented in the flow chart of Fig. 8 . Examples that illustrate the outcome of these processing steps on a vessel that exhibits both phase wrapping and fringe washout can be seen in Fig. 7. The importance for reconstructing the correct pulsatile flow behavior is shown in Fig. 9 . Note that only missing data has been reconstructed via interpolation, leaving original Doppler values untouched.
In order to study the repeatability of the procedure, we evaluated three heart cycles of five vessels from three different scans, recorded at two different days with a healthy individual. The pulse rate on both days was similar (72 and 74). The recording frequency of the circular scans that consisted of 3000 A-scans was 10 Hz at an acquisition speed of 30kHz. All images were processed according to the above mentioned pipeline and the flow measurements were corrected for the angles.
The systolic and diastolic velocity values are extracted from time series as shown in
Fig. 9(b). The results that have been achieved on five selected vessels (see Fig. 3) are summarized in Table 1 . In addition, the coefficient of variation c (standard deviation/mean) is given. “A” stands for artery, “V” stands for vein:
With a recording frequency of 10 Hz of the circular scans, we are at the limit of sufficiently resolving the systole and diastole phases of the heart. This influences the exactness of extracting the correct minima and maxima and therefore directly affects the precision of calculating RI and PI. This is the reason why the results for PI and RI are less reproducible than those for the average velocity which are less affected by the correct minimum and maximum sampling during a heart cycle (see Fig. 9).
The finite sampling of the volume might reduce the velocity precision if the circular scan is located just between two positions of the reference volume where the angle is known. Small errors in the registration might then lead to a wrong assignment of angles. A possibility to cope with this problem might be to record two volumes with orthogonal fast scan axes and register them to each other. After interpolation this might lead to a better sampling of the vessels, allowing for the extraction of more precise angles.
The sampling of the volume also has a large effect on the accuracy of the angles in regions of great variation of the vessel orientation with respect to the optical axis, e.g. near a turning point. This is especially critical as small variations of angles close to 90 degrees lead to large variations of 1/cos α. Small errors in the angles, therefore, greatly affect the corrected flow values. For the vessels that have been examined in this study, the angles were between 78 and 83 degrees.
In general the recorded reference volume should exhibit as much useful phase information as possible, as the manual segmentation of the vessels is done on the Doppler data. Useful phase information shows clear Doppler signature above bulk noise without being cancelled by fringe washout. In order to keep the recording time for a full volume reasonable, the recording speed and volume sampling density need to be well adjusted. Large recording speed results in lower SNR and lower velocity sensitivity keeping the same sampling density. Increasing the sampling density again increases the recording time. A compromise has to be found between recording speed and the tolerance of motion artifacts.
It might happen that the orientation of a circular scan does not fit the orientation of the reference volume if the position of the eye with respect to the incident beam changed between the recording of the reference volume and the circular scan. However, this would immediately be visible during the registration: most clearly the surface shape of the circular scan would not match that of any virtual scan extracted from the reference volume. In general it would be impossible to properly register the scan to the volume. The similarity of the surfaces as given in Fig. 4(c) after the registration is a good indicator for the correct orientation of the circular scan with respect to the reference volume.
Furthermore, registration errors can occur if the actual scan pattern is distorted beyond a simple elliptical shape. An example of this can be seen in Fig. 10 . As can be seen, the intensity profiles of the circular and virtual scan drift apart at the end of the scan. The fact that two-thirds of the profile are aligned correctly leads to the assumption that the scan might have an egg shape. This assumption is also supported by the fact that different optimization runs match either side of the scan perfectly. An egg shape basically consists of two different ellipses joined at the equator with the same axis of rotational symmetry. However, such shape was rather an exception during the measurements, which is why we did not include it in the optimization routine.
We introduced a motion-artifact stable algorithm for extraction of the correct Doppler velocity within retinal vessels. The method was based on referencing circumpapillary scans to a 3D volume by 2D/3D registration. The correct local vessel orientation within the 3D reference volume was segmented once and the extracted angles could then be accessed via a look-up table. The method is capable of resolving pulsatile flow over time for individual retinal vessels. In addition, we corrected single phase-wrapping of the flow and fringe wash-out by interpolation based on the assumption of a parabolic flow profile. The algorithm was successfully tested on vessels of a healthy individual. The repeatability of the extracted velocity was tested showing a variation of only 3-8%. We further tested the reliability of the resistance index and pulsatility index extraction. The indices varied by 8-21%, which is due to our finite temporal sampling. We believe that the present work is an important step towards reliable perfusion assessment using Doppler OCT. The possibility to have reproducible quantities for characterizing retinal perfusion will play an important role in future retinal diagnosis.
We acknowledge financial support from the European Commission Seventh Framework Programme (FP 7) HEALTH program (grant 201880, FUN OCT).