1. Introduction

Fourier domain optical coherence tomography (FDOCT) has nowadays become the method of choice for clinical ophthalmic imaging of the human retina [1, 2]. The reasons are its outstanding sensitivity for high resolution imaging as well as achievable imaging speeds of typically several 10,000 depth profiles/second [3–5]. Recent advances in swept source technology brought an increase in imaging speed of a factor of 10[6]. Currently such sources operate in the near infrared region above 1µm, with limited sweeping range. The spectral range is still larger for spectrometer-based systems, where a depth resolution of 3µm has been demonstrated [7, 8]. The introduction of sensitive CMOS technology brought a similar improvement in imaging speed comparable to swept-source systems. Retinal images at depth profile rates of more than 300,000 scans/second have been presented at reasonable sensitivity [9]. The increased imaging speed is especially interesting for functional extensions of FDOCT, in particular for optical Doppler tomography (D-FDOCT) [10–13]. Perfusion is an important parameter and indicator of tissue health. Several studies pointed to the sensitivity of perfusion figures to pathologic retinal changes already at an early state, such as [14–16]. The sensitivity and speed advantage of D-FDOCT offered for the first time the possibility to image the 3D perfusion structure of the retina [17]. Completely non-invasive optical retinal angiography came into reach [18]. Different approaches have been developed to extract tissue perfusion from surrounding tissue. Lately, techniques based on B-scan analysis of Doppler data achieved separation of flow components from static tissue showing impressive results of vascularization fine structure [19–21]. Nevertheless, their performance strongly depends on the lateral spatial spectrum and on the oversampling ratio.

All methods so far suffer from a lack of information about heart beat during data recording. Since a full 3D dataset of adequate sampling density requires several seconds to acquire, a volume stretches over multiple heart cycles and different tomograms correspond to different pulse phases. Hence, it is difficult to get a coherent impression of flow dynamics within tissue volumes. In principle the periodic characteristic of the pulse can be used to resample several volumes, in order to achieve individual 3D sets corresponding to a single heart beat cycle. However, motion artifacts, which are unavoidable for human retinal imaging, leave such task highly challenging. Currently a gated 4D Doppler approach has been only demonstrated in animal studies [22].

The ultra-high-speed FDOCT approach offers completely new perspectives for 4D imaging. In the following we would like to study the performance of this system for D-FDOCT imaging at different speeds, and compare the results to theoretical expectations. Further on, we demonstrate to our knowledge first 4D imaging of human retinal vessels, acquired at a rate of 13.3 volumes per second.

2. Experimental Setup

 

Fig. 1. Optical Setup.

Download Full Size | PDF

The optical setup is displayed in Fig. 1. We employ a superluminescence diode (SLD) of 45nm spectral width centered at 830nm (EXALOS Inc.) resulting in an axial resolution in air of 6.7µm. The light is fiber coupled into a bulk-optics Michelson interferometer. In the sample arm, the light passes a X-Y galvo-scanner unit, and is directed to the eye through a telescope (L1-L2). We used two telescope settings: L1=60mm, L2=40mm for standard imaging, and exchanging L2 and L1 for high-resolution imaging. The respective beam diameter at the cornea is 2mm and 7mm. For all experiments the power of light incident on the eye was measured to be 550µW, which is in conformity with the IEC laser safety standards [23]. The reference arm is equipped with a dispersion compensating prism-pair and water chamber (25mm), focusing optics, and a mirror mounted on a linear translation stage. The light from sample and reference path is combined at the exit of the interferometer and guided through a single mode fiber to the spectrometer. The latter consists of a volume diffraction grating (Wasatch, 1200lines/mm), an achromatic camera lens (f=150mm), and the high speed CMOS line sensor (Basler Sprint spL4096-140km), with 4k pixels. We image the spectral half width of the SLD of 45nm onto 814 pixels of the sensor resulting in a full depth range in air of 3.1mm. The effective maximum speed of the detector depends on the chosen number of horizontal pixels. In order to compare the performance at various imaging speeds, as well as for 4D imaging, we chose 1152 pixels allowing for a maximum line rate of 200,000 Hz. It should however be mentioned that this setting apodizes already the spectrum of the SLD.

 

Fig. 2. (a) spectral intensity of SLD. The yellow curves show the used apodized spectrum. (b) Coherence function for full spectrum (dotted line) and for the used spectrum (solid line). The intensity is normalized to the maximum of the original coherence function (dotted line). (c) Sensitivity decay with depth (red line 20kHz, black line 200kHz). (d) Sensitivity decay with line rate (black), and actual exposure time (blue).

Download Full Size | PDF

As seen from Fig. 2(a–b) the resulting degradation of depth resolution as well as SNR is negligible. The sensitivity was measured using a mirror as sample and a neutral density filter with optical density 3.0 in front of the mirror resulting in a SNR offset of 60 dB. The noise level was obtained from the standard deviation of the noise floor σ_N close to the mirror signal amplitude SM. The measured sensitivity Σ is then Σ(dB)=60dB+SNR(dB), with SNR calculated as SNR=(S_M-N)²/σ² _N, and N being the mean value of the noise floor. The sensitivity decay is displayed in Fig. 2(c) for two camera line-rate settings, 20kHz and 200kHz. It is in both cases -9dB over 2mm, hence there is no observable pixel cross talk change with line rate. Figure 2(d) displays the measured sensitivities for used camera line rates. Since sensitivity scales with exposure time rather than line rate, we plotted the actual exposure time, which does not depend linearly on line rate. There is -11dB change in sensitivity from 20kHz to 200kHz as should be expected, since the exposure times change by an order of magnitude from 48.3µs to 3.3µs. The sensitivity might be partially recovered using alternative interferometer configurations employing fiber circulators or different splitting ratios. [9]

3. Theory

3.1. Doppler Data Extraction

FDOCT signal reconstruction gives easy access to interferometric phase information as a function of depth. This information has per se no absolute value, but the difference yields an absolute value for small axial structural changes in the nanometer range. For unambiguous velocity determination the phase difference Δϕ needs to be confined to [-π, π]. The unwrapped phase difference ΔΦ will be Δϕ for |Δϕ|≤π and Δϕ-sign (Δϕ)2π for |Δϕ|>π. Since the time difference between successive A-scans is well defined by the detector line rate, one can immediately assign a velocity to the unwrapped phase difference ΔΦ [24]

where λ is the central light source wavelength, T the A-scan period, m the number of averaged lateral A-scan phase differences, and α the angle between velocity and observation direction. There are different relations in the literature to extract the velocity. We chose the above form (Eq. 1), since it yields reliable values for large velocities close to the limits. Equally, there have been different approaches to extract the flow angle α [13, 25–26]. Probably the most precise methods increase the number of illumination directions—ideally sampled synchronously—in order to solve for the angle ambiguity [27–29]. An alternative approach is to define angle independent perfusion parameters.

As mentioned above the maximum flow that can be unambiguously quantified corresponds to a phase difference of ±π, i.e., v_max=±λ/(4T). The minimum speed is given by the phase noise ΔΦ_err present in the system, i.e., v_min=λΔΦ_err/(4πT). A detailed discussion of factors that define the velocity bandwidth is presented in section 3.3. The velocity dynamics can be written as 2v_max/v_min=2ΔΦerr/π.

3.2. Background Velocity Correction

Any scanning tomography method suffers from motion artifacts. Even if the time for one scanning direction is as fast as a few milliseconds, highly phase sensitive detection will still pick up phase offsets due to back ground motion. Figure 3(a) shows a typical retinal Doppler tomogram taken at 20kHz with well visible background motion. Different approaches have been suggested for correction of those artifacts. We adapted the method of Makita et al. [18] with a slight modification to stabilize the correction also in the presence of large retinal vessels. We will not revisit the original method in detail, since it is well explained in the literature. Roughly, it corrects the background motion by calculating the shift of the Doppler histogram maximum, which should represent the Doppler offset of the static structure. The correction is performed within each window with its center sliding across the B-scan. If the window size is too small compared to the vessel size, one picks up artifacts, since the main contribution to the histogram maximum originates from the flow signal (see arrow in Fig. 3(b)). If the window is large the motion correction according to [18] will be increasingly time-consuming. Our modification of the algorithm is based on the assumption that background motion gives rise to a continuous change of the histogram maximum. In this case we can fit a polynomial function through the correction shift values, and in a second step detect artifacts as large deviations from this curve (Fig. 3(c)). They are corrected by cropping the phase difference values to within a maximum deviation band from the fitted polynomial (dashed lines in Fig. 3(c)). The order of the polynomial function is chosen to be of low order to avoid that the Doppler shift curve follows sudden jumps. A polynomial of order three gave stable results with respect to artifact free background motion correction. The final result is shown in Fig. 3(d). The control parameters of the modified algorithm are the maximum allowed deviation from the polynomial as well as the window size. Both are chosen empirically and held fixed for the full 3D datasets.

 

Fig. 3. (a) Phase difference tomogram at 20kHz covering 4deg near the ONH. (b) Motion correction according to Makita et. al. [18] with small window size. Arrow points to an artifact. (c) Modified motion correction by limiting phase difference deviations (blue line) from fitted polynomial (red line) to area within dashed lines. (d) Result after modified motion correction algorithm.

Download Full Size | PDF

3.3. True D-FDOCT velocity dynamics

Imaging flow with spectrometer based FDOCT has an important drawback: the fringe signal corresponding to moving structure will shift during camera exposure resulting in fringe washout and degradation of signal-to-noise ratio (SNR). With increasing speed relative to line rate the SNR attenuation follows approximately a sinc-function (Fig. 4(b)) [17, 30]. It should be mentioned that for larger velocities the attenuation curve deviates from the sinc-function due to large spectral bandwidths employed in OCT-systems that cause the attenuation curve to never completely drop to zero. In resonant Doppler imaging the attenuation curve maximum can be shifted via reference phase shifting to enhance moving structures and suppress static tissue components [17, 25]. If we increase the line rate and shorten the exposure time the attenuation curve will be broadened (Fig. 4(b)) and a larger velocity range stays above the SNR limit of the system (Note that the curve is a relative curve with its maximum representing SNR if the associated moving structure was at rest). Let us now discuss in more detail the limiting factors for the actual velocity range. At the lower end the minimal resolvable Doppler velocity is determined by the system phase error. Usually the phase error is calculated from the standard deviation of temporal phase fluctuations at a fixed lateral position on a sample mirror. Setting up this configuration we measured a phase error of 0.032 rad. This error stayed the same for different camera line rate settings, keeping the same SNR. It has been shown that such phase error is not representative for phase fluctuations present, if scanning across a scattering sample is employed. Following [31] the phase error is independent of camera line rate and related to the inverse lateral oversampling factor g=w/(sN), where w is the geometric tomogram width, N the number of A-scans per tomogram, and s the spot size. In our standard experimental setting (cf. section 2 and 4) for measuring the optic nerve head region the width of the full transverse tomogram was 2mm with a spot size of 20µm consisting of 1000 A-scans. The resulting sampling factor of 0.1 gives rise to a theoretic phase error of 0.4 rad (cf. Fig. 6 in [31]).

 

Fig. 4. (a) Histogram of phase differences for static structure at various line rates. The widths are almost identical. (b) Attenuation curves for different line rates depending on velocity. Dotted lines indicate the unambiguous maximum velocity limit (vertical line) and visibility limit (horizontal line).

Download Full Size | PDF

We suggest estimating the experimental phase fluctuations ΔΦ_err from the width of the phase difference distribution histogram of a lateral tomogram section that contains mainly static structure. Figure 4(a) displays the phase difference histograms for different line rates. The theoretic value for ΔΦ_err is well reproduced by half of the average full-width-half-maximum of the histogram plots, which was 0.41 rad. Comparing this value to the original static phase error of 0.032 rad we observe a discrepancy by an order of magnitude which has a direct impact on the minimum resolvable flow speed. We therefore recommend for future D-OCT systems to use the suggested performance figure for a better comparison across different D-OCT platforms.

At the high velocity end we have two limiting factors. Firstly, unambiguous flow quantification sets a limit of ± π for the maximum phase difference (cf. section 3.1). Secondly, the visibility of structure in motion is determined by the maximum achievable attenuation, which depends on the assumed SNR of the structure, if it was at rest [17]. The latter depends on various factors such as position of the structure within the system depth range, position of structure within the tissue, depth of the structure within a vessel due to strongly scattering red blood cells, and last but not least on the exposure time. Since the SNR issue is obviously difficult to generalize, we suggest using the average SNR of static structure as reference. Those values are given in Table 1 for different detector line rates. If we plot those values together with the attenuation curves (dotted horizontal lines in Fig. 4(b)) we obtain the largest velocity of structure in motion given the associated average SNR at that line rate. From the discussion above follows that this is not a sharp limit since it depends on all factors that impact on SNR. Nevertheless, we observe, that in average the limit of unambiguous flow calculation (vertical dotted lines in Fig. 4(b)) is stronger than the SNR limit. SNR might ultimately also reduce the velocity sensitivity, since the system phase error can be expressed as ΔΦ=1/√SNR [30, 32]. Hence, if we had theoretically increased the oversampling factor, we would eventually hit the SNR boundary. In practice the limit set by phase errors introduced through scanning are dominating. The histogram plot in Fig. 4(a) verifies this fact: if SNR would be dominating then the phase error (histogram width) would increase with decreasing exposure time.

Table 1 summarizes the performance figures of our setup for different detector line rates. We distinguish between the high stability case where no scanning is employed and a mirror was used as sample (wo/scan) with ΔΦ_err=0.032 rad to the scanning case (w/scan) with ΔΦ_err=0.41 rad for our experimental configuration. For the highest velocity we gave the maximum velocity figure for the visibility of structure in motion (SNR) and the value for the highest non-ambiguously quantifiable velocity (na). In principle phase unwrapping would allow to exploit the full velocity range defined by the visibility of structure (SNR). The last column gives the inverse of the quantifiable velocity bandwidth (see section 3.1).

Table 1. Standard D-FDOCT performance figures for different array detector line rates

View Table

The above table refers to the case of keeping lateral oversampling constant. One could alternatively increase the lateral oversampling resulting in smaller phase noise. Although one looses then the advantage of increased imaging speed one might gain higher dynamic velocity range. The loss in velocity sensitivity is a major obstacle for quantifying flow at high acquisition speeds. At 100 kHz we exceed already typical measurable velocities present in retinal vessels due to their orientation almost perpendicular to the detection axis. From these considerations follows that higher imaging speed will not in general enhance quantitative Doppler tomography. Nevertheless, if the speed is sufficiently high, such that volume time series can be taken to sample the full heart cycle, higher speed opens a new dimension in perfusion tomography. In combination with high resolution OCT we will show how to assess the full 4D dynamics of the retinal capillary bed.

4. Results and Discussion

4.1 Performance of 3D D-FDOCT at different acquisition speeds

To evaluate the performance of 3D D-FDOCT at different acquisition speeds we recorded volumes of 1000 A-scans×200 B-scans with line rates of 20kHz, 60kHz and 100kHz covering a range of 4deg×4deg across the ONH (Fig. 5, (View 1), (View 2), (View 3)). For these experiments we used a beam diameter at the cornea of 2mm with a spot size of 20µm on the retina and a sampling factor of 0.1. We chose to image the ONH in our experiment, because it yields the larges dynamic range of blood flow, due to the varying orientations of the vessels in this region. An additional volume of 1000 A-scans x 200 B-scans covering a range of 1.2deg×1.2deg at the fovea was recorded at a line rate of 160kHz to try to quantitatively resolve flow within small retinal vessels (Fig. 6, (View 4)). For this experiment we enlarged the beam diameter at the cornea to 7mm to improve transversal resolution, resulting in a spot size of 5µm on the retina and a sampling factor of 0.12. Axial movement artifacts within the volumes were corrected using correlation based image registration. The phase difference tomograms were additionally corrected for background movement as mentioned above.

Figure 5 shows for each of the imaging speeds a volume with corresponding cross-sections at significant positions. The shown 3D datasets were saved in the portable network graphics (png) file format. Besides the standard RGB channels this format contains an additional transparency channel (α-channel). Our files contain the phase difference or Doppler tomograms in the RGB channels, with green for static structure and red and blue corresponding to positive and negative motion. The transparency. channels of the png-encoded images were used to store the intensity tomograms. Tuning the transparency (α-channel) allows then to suppress the noise of the phase tomograms and to enhance structural details. Fig. 5J-L illustrates the option to choose a different color map within OSA ISP in order to improve contrast of Doppler information. Alternative color maps are available and can be uploaded during an OSA ISP session.

 

Fig. 5. Performance of 3D D-FDOCT at different acquisition speeds covering a patch of 4deg×4deg across the ONH. (A, (View 1), D, (View 2), G, (View 3), J) Volumes saved in the png file format containing D-FDOCT tomograms in the RGB channel and intensity tomograms in the α-channel. The volumes were recorded at line rates of 20kHz(A and J), 60kHz(D) and 100kHz(G). (B, E, H, K) corresponding en-face cross-section at a significant position. (C, F, I, L) corresponding fast axis cross-sections (B-scan) at significant positions. (J, K, L) demonstrate the possibility of a different color map on a 3D dataset acquired at a line rate of 20kHz. Dashed arrows in (B, E, H, K) show vessels with high flow speeds and large inclination. Solid arrows in (B, E, H, K) show vessels with lower flow speeds and almost perpendicular orientation to the detection axis. Interferometric data of measurements in Fig.5A, D, G are provided as download.

Download Full Size | PDF

For 20kHz (Fig. 5(A-C), (View 1)) we receive good results for smaller vessels orientated almost perpendicular to the detection axis (Fig. 5B, solid arrow). Vessels coming out of the ONH more steeply appear hollow due to fringe blurring (Fig. 5B, dashed arrow). The volumes recorded at a line rate of 60kHz show the best results concerning D-FDOCT signal visibility (Fig. 5(D-F), (View 2)). The dynamic range of detectable flow speeds at this line rate matches the flow speeds within the volunteer’s ONH region very well. Vessels lying rather perpendicular to the detection axis (Fig. 5E, solid arrow) as well as vessels more parallel to the detection axis (Fig. 5E, dashed arrow) both yield D-FDOCT information. The volumes acquired at 100kHz (Fig. 5(G-I), (View 3)) show anticipated results. All vessels appear full; however only within the largest vessels and for vessels almost parallel to the detection axis (Fig. 5H, dashed arrow) the D-FDOCT analysis delivers quantitative results. For the smaller vessels and vessels lying perpendicular to the detection axis (Fig. 5H, solid arrow), the time between two adjacent A-scans is too short and the velocity resolution too small to observe a phase difference.

We experienced the same predicament with the measurements in the fovea region. Even though imaging at such high speed reduces motion artifacts to a minimum and therefore enables imaging of very small details as the capillary bed, we are not able to resolve the slow capillary flow. Fig. 6 shows intensity en-face views and D-FDOCT en-face views of vessels in the ganglion cell layer (GCL) (Fig. 6(A, B)) and capillaries in the inner nuclear layer (INL) (Fig. 6(C, D)). These views were extracted from a volume of 1000 A-scans × 200 B-scans × 576 axial pixels (View 4), acquired at a line rate of 160kHz, resulting in an acquisition time of 1.25 sec.

 

Fig. 6. Flow within small retinal vessels. Intensity en-face view of small vessels within the GCL (A) and the capillary network within the INL (C) (View 4). (B and D) corresponding D-FDOCT en-face views covering a patch of 1.2deg × 1.2deg across the fovea.

Download Full Size | PDF

4.2 4D Micro-FDOCT

Large imaging speed opens completely new perspectives for imaging small retinal details. Fast acquisition allows capturing fine retinal details by strongly reducing transverse motion blurring. For the following experiments we used the enlarged beam diameter at the cornea of 7mm, resulting in a spot size of 5µm on the retina. We operated the camera at 200.000 kHz Fig. 7(A) demonstrates the resolution capabilities visualizing retinal cone mosaic without the application of adaptive optics[9, 33]. The volume of 500 A-scans × 200 B-scans covering a range of 0.4deg × 0.4deg is recorded at a foveal eccentricity of 4deg. As we have seen in the previous section high speed standard D-FDOCT will not be able to quantitatively resolve flow within small retinal vessels, but it will help to image with good resolution and small motion artifacts the retinal capillary bed (Fig. 6(A, C), (View 4)). In order to observe dynamic perfusion changes we reduce the lateral sampling to 150 A-scans × 100 B-scans, resulting in a volume rate of 13.3 per second. This rate is sufficient to resolve the heart cycle. We recorded a retinal patch of 0.6deg × 0.4deg close to the fovea at an approximate eccentricity of 1deg. Axial movement artifacts within the volumes were corrected using correlation based image registration, transversal movement between adjacent volumes was corrected manually. The 4D movie (Media 1) shows a pass-through through the retinal layers, the additional temporal information enables direct visualization of flow within the capillary network. The video starts out with an en-face view onto the nerve fiber layer (NFL) then passes through the NFL to reveal a vessel in the GCL (Fig. 7(B), (Media 1)). First the flow is poorly visible from outside the vessel, but when another few slices were taken off one can clearly see the flow within the vessel. It then continues to cut off further slices to uncover first the capillary network superior then anterior to the INL. Even within some of these capillaries flow is noticeable. After passing through the photoreceptor layers the choroid becomes visible. At the end the complete volume is rebuild to show the patch with all layers.

To the best of our knowledge it is the first time that retinal flow has been visualized within a full 3D volume over time.

 

Fig. 7 4D Micro-FDOCT.(A) cone photoreceptors imaged without the application of adaptive optics or active eye tracking. Size: 0.4deg × 0.4deg at 4deg eccentricity. (B) volume with removed NFL covering 0.6deg × 0.4 deg at 4 deg eccentricity showing a vessel in the GCL (Media 1).

Download Full Size | PDF

7. Conclusion

In conclusion we demonstrated the performance and limitations of ultra-high-speed spectrometer based FDOCT for resolving retinal perfusion. We point out relevant figures to characterize velocity performance across different D-OCT platforms. Standard Doppler evaluation will already fail to quantify retinal perfusion for most of retinal vessels at line rates above 100kHz. Nevertheless line rates above 200kHz allow recording volume time series that reveal perfusion dynamics on the intensity level. Although flow quantification has not been demonstrated within such time series, it offers for the first time the observation of pulsation phase-coherent volumetric perfusion images. Such information will certainly impact on the way retinal flow beds are understood in their interaction with different retinal structures and layers and might eventually lead to a better understanding of retinal function and patho-physiology.