1. Introduction

Spectral domain optical coherence tomography (SDOCT) has demonstrated clinical potential for in vivo high-resolution and high-speed imaging of biological structures [1, 2]. Advances in Doppler SDOCT (DOCT) have demonstrated several image acquisition schemes that enable real-time, high-resolution, volumetric display of blood flow maps [3–8]. Recently, several methods have been introduced which focus on utilizing spatial frequency modulations across lateral scans to resolve moving scatterers within a B-scan [9–12]. Techniques have also been presented which use reference delay modulations to specify a velocity measurement range [13]; or uses temporal frequency shifts to improve velocity resolution in low signal conditions [14]. The theoretical bases behind these techniques are closely related to those used in complex conjugate resolved imaging [15–23]. A subset of these techniques [9, 10, 12] have demonstrated improved flow detection over conventional DOCT, resolving in vivo microvascular networks. Compared with DOCT, however, this subset of techniques has sacrificed velocity resolution for improved flow sensitivity, providing only flow direction without speed.

Spatial frequency flow detection techniques [9, 10, 12] can be considered optical analogs to power Doppler (PD) ultrasonography [24–26]. Developed as a method of improving the sensitivity of Doppler ultrasound, the analog to DOCT, PD reports the power of the Doppler signal within specified frequency windows instead of the mean frequency shift. The advantage arises from the representation of the power spectrum of random phase noise. Since the noise in the power spectrum is uniformly low, random phase variations can be filtered out by raising the sensitivity threshold above the noise floor. The Doppler signal in PD is represented as an integral of the power spectrum, which improves the sensitivity and detection range of moving scatterers at the expense of eliminating velocity information. PD is relatively insensitive to Doppler angle and phase wrapping since these factors only modify the distribution of the Doppler power spectrum, but the total integrated power remains constant. The Doppler signal in PD is separated from non-moving components by filtering out all power spectrum components without Doppler shifts, thus only imaging moving scatterers. The resulting PD signal is related to the number of moving scatterers producing the corresponding Doppler shifts. Algorithmically, PD is equivalent to the single-pass flow imaging SDOCT (SPFI-SDOCT) technique we recently demonstrated [12].

Spatial frequency filtering techniques [9, 10, 12], while able to provide 3D flow maps with improved acquisition speed and sensitivity compared with DOCT, lack velocity-resolved blood flow information provided by techniques such as DOCT [4] and laser Doppler velocimetry (LDV) [27]. In ophthalmology, noninvasive quantification of blood circulation in tissues can facilitate the description of retinal vascular changes prior to and during ocular and systemic disease. Here, we demonstrate an improvement to SPFI-SDOCT which allows for velocity-resolved flow mapping using a similar modified Hilbert transform algorithm. The spatial frequencies of moving scatterers detected by SPFI are windowed such that only a small spatial frequency range, corresponding to a finite velocity range, is reconstructed. The window is then scanned across all spatial frequencies, forming a stack of B-scans each corresponding to a different velocity. The stack can then be summed across the entire velocity range to reconstruct a single velocity-resolved B-scan. OCT volume datasets can be processed to yield three-dimensional, velocity-resolved flow maps featuring the flow detection sensitivity improvement of SPFI, and the velocity resolution of conventional DOCT. Flow velocities calculated using this technique is validated against measurements taken with DOCT in flow phantoms. Velocity-resolved SPFI is then used to reconstruct flow maps of in vivo human retina estimating, for the first time to our knowledge, total foveal blood flow using OCT.

2. Theory

The theoretical basis for SPFI-SDOCT [12], and similar techniques [9, 10, 13, 17], have been developed previously. Briefly, for a recorded interferometric B-scan, I[k, x], the spatial Fourier transform along the lateral scan dimension (Fig. 1(a)) yields the spatial frequency content

Here, both moving, V _±[k, u+f _{D, ±}], and non-moving, R′[u, k], scatterers in the sample are imaged (Fig. 1(b)), where f _{D, ±} represents the Doppler frequency shift associated with the axial components of scatterer motion. All scatterer motion and associated flow refer to the axial components of motion.

In SPFI-SDOCT, we recognize that the spatial frequencies of moving and non-moving scatterers do not overlap at spatial frequencies above the non-moving scatterer bandwidth (Fig. 1(b)). An analytic signal for the spectral interferogram can be obtained by applying a Heaviside function (H[u-f_T]), frequency-shifted outside of the structural bandwidth (Fig. 1(c)), and then inverse Fourier transforming the result (Fig. 1(d)). Application of this modified Hilbert transform (HT*) enables bidirectional flow imaging by windowing Eq. (1) to yield

where α represents the fractional portion of bidirectional flow with Doppler frequencies outside of the spatial frequency bandwidth of non-moving scatterers (Fig. 1(c)). This threshold frequency (f_T) defines the minimum detectable velocity in SPFI-SDOCT and is related to the spatial correlation of sequential A-scans [12]. The spatial oversampling and velocity resolution can be related by combining the velocity-related Doppler frequency shift with the spatial frequency resolution, thus yielding a velocity resolution

where λ ₀ is the center wavelength, τ is the integration time, and θ_D is the Doppler angle between the scanning beam and the direction of scatterer motion. n is the index of refraction, L is the lateral scan length, D is the number of A-scans acquired across the lateral scan, and w_o is the scanning beam spot size. Eq. (3) shows that velocity resolution in SPFI increases and the maximum detectable velocity decreases linearly with increased spatial oversampling.

Since velocity depends linearly on frequency, application of a spatial frequency window will necessarily filter out all velocities not traveling at velocities described by

Here, W[u] represents the spatial frequency FWHM bandwidth of the Gaussian window. There is, however, a tradeoff between the spatial frequency window width and the resulting spatial resolution of the velocity-resolved intensity image. For a Gaussian window in frequency-space, the resulting spatially resolved vessel map will be convolved with a Gaussian in the B-scan dimension with FWHM bandwidth W[x]=4ln2/(πW[u]). Thus, to avoid loss of lateral resolution in the velocity-resolved vessel maps, the spatial frequency window bandwidth, W[u], needs to be constrained such that the associated lateral blurring function, W[x], does not exceed the scanning beam spot size, w_o. The velocity window limit can be represented as the bandwidth W[u]≥4ln2/(πw_o). Given the discrete sampling parameters of SPFI-SDOCT, the minimum W[u] before loss of resolution is a factor of 8ln2/π greater than the spatial frequency sampling rate. Combining this result with Eq. (3) gives the minimum resolvable velocity resolution as

 

Fig. 1. Flow chart of velocity-resolved SPFI-SDOCT signal processing. (a) Lateral Fourier transform of raw SDOCT spectral interferogram B-scan yields (b) spatial frequency of stationary scatterers centered around DC, and spatial frequency of moving scatterers shifted by their respective Doppler frequencies. (c) Applying a frequency-shifted Heaviside step function, spatial frequency windowing, and inverse Fourier transforming each frequency range recreates (d) the analytic interferometric signal. (e) Spectral inverse Fourier transform of the analytic interferometric signal maps depth-solved reflectivities of moving scatterers for each corresponding velocity range into a datacube for each B-scan. Bidirectional flow is mapped onto opposite image half-planes. Summing the datacube across all velocity ranges creates velocity- and depth-resolved B-scans.

Download Full Size | PDF

Finally, a datacube, I[x, z, v], for each B-scan can be formed by shifting the spatial frequency window across the bandwidth of moving scatterers (Fig. 1(c)), inverse transforming the result back to functions of I[k, x] (Fig. 1(d)), and then performing a spectral inverse Fourier transform (Fig. 1(e)). Each v-slice represents a velocity range given by Eq. (4), and the datacube can be summed across the v-dimension to create a single velocity-resolved B-scan. Similar to conventional DOCT, velocity wrapping occurs as the spatial frequency content of moving scatterers wraps across the Nyquist sampling upper limit. In this case, the scatterers are mapped to the opposite image half-plane and are therefore represented as moving in the opposite direction. Phase unwrapping techniques, similar to those used for DOCT [28], can be applied for a singly wrapped velocity profile. Higher spatial frequencies that wrap across Nyquist and have components greater than the threshold frequency, f_T, will necessarily lose their velocity components due to SPFI windowing. This fundamentally limits the resolvable spatial frequency range to twice the Nyquist frequency.

3. Methods

Velocity-resolved SPFI-SDOCT was implemented on a high-speed SDOCT retinal imaging system (Fig. 2) employing a light source with a central wavelength of 859nm and a FWHM bandwidth of 99nm (Superlum, Ltd.). The sample arm was a modified slit lamp equipped with scanning galvanometers and relay optics for retinal imaging of subjects. The retinal scanner optics were designed for a 15–20µm transverse resolution, as limited by the optics of the eye, across a 12×12mm field. The reference arm was dispersion compensated using a water cell and matched optics, and interferometric signals were captured using a 2048 pixel line-scan camera (e2v, Ltd.). Custom software (Bioptigen, Inc.) performed data acquisition, archiving, and real-time processing and display of image magnitude. Using a 700µW sample beam, the SNR measured near DC was 110dB with an axial resolution of 4.72µm in tissue and a 6dB falloff at 0.8mm. DC removal, k-space resampling, and flow imaging using the modified Hilbert transform algorithm [12] were computed during post-processing using Matlab (MathWorks, Inc.). Vessels and structure were visualized using Amira (Visage Imaging, Inc.) and OSA ISP (Kitware, Inc.).

 

Fig. 2. Velocity-resolved SPFI-SDOCT retinal system. The sample arm is a modified slit lamp equipped with scanning galvanometers and relay optics for convenient patient retinal imaging. SPFI windowing and velocity-map reconstruction were calculated in post-processing

Download Full Size | PDF

Velocity-resolved bidirectional flow imaging was validated and compared with conventional DOCT on a flow phantom. Two glass micro-capillary tubes (1.5mm outer diameter, 0.6mm inner diameter) were connected using silastic tubing to a syringe pump (Harvard Apparatus) and pumped with 1% liposyn at 10µL/min, 20µL/min, 30µL/min, and 40µL/min. The micro-capillaries were then positioned adjacent to each other on an angled stage such that fluid in the tubes flowed in opposite directions in a B-scan cross-section (Fig. 3(a)), simulating bidirectional flow. B-scans of the phantom were acquired across a 2mm scan range with 2500 A-scans/frame for SPFI and 1000 A-scans/frame with 4 sequential A-scans at each lateral position for DOCT. The SPFI dataset was laterally oversampled compared to DOCT because velocity resolution increases as a function of spot-size overlap on the sample (Eq. (3)). Both datasets were acquired with an integration time of 50µs. At these sampling parameters, the total imaging time for SPFI is a factor of 1.6 times faster than that of DOCT. SPFI parameters were chosen to demonstrate velocity resolution at the lower limit of the detection range for a given integration time as a comparison with DOCT. Since the velocity resolution of DOCT is limited by the phase noise of the system instead of spatial sampling parameters, as is the case with SPFI, an appropriate lateral spacing used to minimize scanner jitter. The number of sequential A-scans used is indicative of common DOCT sampling parameters.

 

Fig. 3. Validation of SPFI measured flow velocities. (a) Two micro-capillaries were connected and oriented such that 1% liposyn flowed in opposite directions in a B-scan cross-section. (b),(c) Velocity profiles measured using DOCT and fit to laminar flow curves. Data was acquired with 1000 A-scans/frame with 4 sequential A-scans for each lateral position with 50µs integration time. (d),(e) Velocity profiles measured using SPFI and fit to laminar flow curves. Data was acquired with 2500 A-scans/frame with 50µs integration time. (f),(g) Flow measured using both Doppler and SPFI were compared for both negative and positive flow, solid line represents theoretical flow rate.

Download Full Size | PDF

A threshold frequency was determined which filtered out all spatial frequencies of non-moving scatterers. A Gaussian window was then moved across the remaining spatial frequencies to velocity-resolve the B-scan for flow rate measurements. The shifted window was set such that the velocity-range resolved had a FWHM and shifted at increments of 24.3µm/s. Given the oversampling parameters and threshold frequency window used, the magnitude of the total detectable velocity range for the axial components of positive and negative velocities was 0.61-11.53mm/s. Velocity-resolved scatterer information was then overlaid onto structural B-scans for visualization (Fig. 3(a)). DOCT volumes were processed using standard phase-difference methods [6]. Velocity profiles for both SPFI and DOCT were fit to laminar flow curves, and measured capillary cross-section and flow rate were calculated for both imaging methods and compared.

In vivo microvessel imaging was demonstrated in normal human retina. First, a 10×10mm OCT volume of the subject was acquired, allowing for the reconstruction of a standard OCT summed voxel projection (SVP) to use as an atlas to locate smaller volumes imaged using SPFI (Fig. (4)). Several 2×2mm volumes were then densely-sampled using SPFI parameters (2500 A-scans/frame, 100 frames/volume, 100µs integration time − 25s total imaging time) at several locations across the macula (Figs. 4(a) (View 1), 4(b) (View 2), and 4(d) (View 4)), including landmarks such as the fovea (Fig. 4(c) (View 3)) and optic nerve (Fig. 4(e) (View 5)). Given the sampling parameters, threshold frequency, and assuming a 20µm scanning beam spot-size on the retina, the detectable velocity range for axially moving scatterers was 0.45-8.64mm/s. Parameters were set for a lower velocity range, compared to flow phantoms, to adequately detect slow flow in small foveal vessels (<30mm/s) [27]. The frequency window was set such that the velocity-range had a FWHM of 18.4µm/s and was shifted at increments of 82.8µm/s. Bulk motion correction using the algorithm described in [10] was implemented prior to SPFI windowing. Velocity-resolved B-scans were used to determine vascular size by fitting velocity profiles to parabolic flow curves and then calculating the zero-velocity crossing positions. Finally, blood flow rates and total retinal flow were calculated for a 2×2mm velocity-resolved volume of human fovea.

SPFI acquisition of the entire macula in a single volume dataset would be impractical. Since velocity resolution and lateral sampling density are coupled in SPFI by Eq. (3), an order of magnitude increase in the lateral sampling density, A-scans/frame, would be required. This requirement limits the SPFI sampling volume size by both the available memory in the acquisition software and total imaging time.

4. Results and Discussion

Velocity-resolved bidirectional flow imaging was demonstrated in a liposyn-pumped flow phantom and compared with measurements taken with conventional DOCT (Fig. 3). Velocity profiles for both positive and negative flow were plotted for both SPFI and DOCT datasets and fit to second-order polynomials, characteristic of laminar flow, to determine the tube diameter and maximum flow velocities (Fig. 3(b)–(d)). The DOCT measured velocity profiles correlated stronger with the parabolic flow fit (mean R²=0.96) than those calculated with SPFI (mean R²=0.81). Since laminar flow velocities across the tube decrease continuously towards zero, all components below the SPFI velocity threshold contribute to velocity noise, which becomes more significant for low flow velocities. The velocities measured using SPFI shows a poor parabolic fit and a correlation of R²<0.7 for the lowest pump speed (Fig. 3(d–e) − dark blue), demonstrating the lower detection limit of the modality. This is a result of the spatial frequency overlap of the slow moving scatterers near the edge of the capillary tubes, which is filtered out in SPFI processing. As scatterer velocity increases, the fit improves and is comparable with that of DOCT (Fig. 3(f–g)). Fig. 3(d) also shows a velocity dip across all pump velocities calculated using SPFI. This is a result of a shadowing artifact as seen in Fig. 3(a), which is more pronounced in SPFI and in DOCT (Fig. 3(b)) since SPFI is an intensity sensitive modality, as opposed to phase sensitive DOCT. For each flow rate, velocity profiles were measured three times for both SPFI and DOCT, and the averaged velocity and standard deviation were calculated (Fig. 3(b)–(g)). Flow rates measured using SPFI and DOCT were plotted as functions of pump flow rates (Fig. 3(f)–(g)), showing strong correlations of both SPFI (mean R²=0.99) and DOCT (mean R²=0.99) with theoretical values.

Depth-resolved vessel maps for each volume were first reconstructed and overlaid onto structural OCT data to distinguish macular vasculature (Fig. 4(a)–(e)). Total blood flow measurements were then calculated for a single 2×2mm volume of the fovea (Fig. 4(c) (View 3)). The structural OCT data shows the foveal pit and the associated SPFI vascular map confirmed its location by resolving a circular avascular zone surrounded by a set of terminal capillaries [29]. First, the resolved vessels in the volume were identified as arteries and veins using the flow directionality information calculated by SPFI (Fig. 5(a)). Vessel orientation and Doppler angle were then measured for all 17 resolved vessels, and velocity profiles were measured at a single point for each vessel (Fig. 5(a) − dots). Doppler angles were measured by calculating the vessel cross-section displacement across sequential B-scans in the 3D datasets [4]. The velocity profiles were then fit to laminar flow curves to determine both the peak velocity and diameter of each vessel. Finally, the vessel size, orientation, and velocity information was used to calculated the total foveal blood flow [4]. A representative B-scan (Fig. 5(b)) taken across the foveal volume (dotted line) is shown with velocity-resolved flow content overlaid on top of the structural image. The velocity profiles and parabolic fits for each of the vessels are included (Fig. 5(c)–(i)) to show the strong correlation between velocities measured using SPFI and their respective laminar flow velocity profiles (mean R²=0.95). Closer inspection of vessel 17 (Fig. 5(i)) shows a blunted parabolic velocity profile, characteristic of red blood cell aggregation in microvasculature, consistent with rheological observations [30–32].

 

Fig. 4. 10×10mm volume of in vivo human retina. (a)–(e) 2×2mm volumes sampled with 1024×2500×100pix at 100µs integration time (25s total imaging time) were acquired at several locations across the macula, (a) (View 1), (b) (View 2), and (d) (View 4), including landmarks such as (c) fovea (View 3) and (e) optic nerve (View 5). Volumes were processed using SPFI, and structure, vessels, and representative B-scans (across dotted lines) with flow overlaid are shown. Detectable vessel diameters ranged from 14µm (fovea) to 120µm (optic nerve).

Download Full Size | PDF

 

Fig. 5. Velocity-resolved 2×2mm volumetric flow map of in vivo human fovea. (a) Fovea vessel map with vessels and flow direction resolved using SPFI. Arrows denote arteries and veins, and dots represent locations where velocity profiles were measured. (b) Representative B-scan depth-slice across the foveal volume (dotted line) with velocity-resolved vessels overlaid onto structural data. Arteries (red) and veins (blue) are labeled for clarity. (c)–(i) Velocity profiles of vessels 11–17, respectively, are fit to laminar flow to determine peak velocities. Individual vessel diameters are measured at the zero-velocity crossings of their respective velocity profiles.

Download Full Size | PDF

A summary of the size, peak velocity, and flow measurements for all vessels identified in the 2×2mm foveal volume are shown in Table 1. The smallest resolvable vessel was 13.64µm, which is at the resolution limit of the retinal SDOCT system. The average measured foveal vessel diameter was ~22µm, and average arterial flow velocity was greater than average venous flow velocity (Table 2), which is supported by similar measurements made using LDV [27]. The detected velocities ranged from 5.97–30.22mm/s, concurrent with human retinal vessels in the corresponding size range, also measured using LDV [27]. Finally, the total arterial and venous flow showed a net inflow of blood into the fovea. This is a result of the presence of unresolved veins in the volume and expected errors from measurements of vessel orientation angle and calculation of diameter, which would significantly impact blood flow calculations.

Table 1. Size, velocity, and flow measurements across a 2×2mm volumetric flow map of in vivo human fovea calculated using SPFI. Arteries and veins (*) are numbered according to those identified in Fig. 5(a). Flow was calculated using vessel sizes and peak velocities measured using the parabolic fits of velocity profiles.

View Table | View all tables in this article

While SPFI allows for velocity-resolved volumetric blood flow imaging with velocity resolution comparable to that of DOCT, there are several advantages and drawbacks to the modality that must be noted. Since velocity resolution and spatial oversampling are coupled, this allows for the flexibility to set sampling parameters to a desired velocity range. This would allow for improved data acquisition speeds over DOCT for imaging of small spatial volumes of moderately high flow velocities, but could potentially require longer scanning times for large scan areas with low flow. As discussed previously, the velocity resolving power of SPFI falls off at the lower limits of velocity detection, as compared to DOCT, but is comparable for high flow velocities. Since velocity-resolved SPFI requires a sliding spatial frequency window, this increases the number of Fourier transforms required by the number of velocity increments desired. However, this increase in computational complexity is essentially trivial since FFT algorithms can be optimized and the transforms can be parallelized. Finally, it has been shown that spatial frequency filtering results in a sensitivity improvement [16], which gives SPFI a detection advantage over DOCT for small vessels at the resolution limit of the imaging system. The overall advantage of SPFI over DOCT is in its improved sensitivity and customizability over velocity resolution, velocity range, and acquisition time in exchange for small increases in computational complexity.

Table 2. Retinal flow statistics in a 2×2mm volumetric flow map of in vivo human fovea calculated using SPFI-SDOCT.

View Table | View all tables in this article

5. Conclusions

We have demonstrated velocity-resolved in vivo volumetric bidirectional flow imaging for the measurements of human retinal blood flow. Velocity-resolved SPFI-SDOCT builds on previously described SPFI imaging techniques, which provide three-dimensional flow imaging without acquiring multiple A-scans at a single lateral position. SPFI measured velocities were validated with velocities measured using conventional DOCT on flow phantoms. Velocity-resolved SPFI was then used for in vivo imaging of several volumes across the macula of healthy human retina. SPFI processing allowed for volumetric visualization of flow without the need for manual segmentation. SPFI also provided flow directionality and velocity, demonstrating the measurement of vessel size, peak velocity, and total foveal blood flow with OCT, consistent with characteristics observed in LDV and rheology.