Volumetric microvascular imaging of human retina using optical coherence tomography with a novel motion contrast technique

: Phase variance-based motion contrast imaging is demonstrated using a spectral domain optical coherence tomography system for the in vivo human retina. This contrast technique spatially identifies locations of motion within the retina primarily associated with vasculature. Histogram-based noise analysis of the motion contrast images was used to reduce the motion noise created by transverse eye motion. En face summation images created from the 3D motion contrast data are presented with segmentation of selected retinal layers to provide non-invasive vascular visualization comparable to currently used invasive angiographic imaging. This motion contrast technique has demonstrated the ability to visualize resolution-limited vasculature independent of vessel orientation and flow velocity.


Introduction
Currently, fluorescein angiography is the gold standard for vascular imaging of the retina. It is commonly used in clinical practice for diagnosis of diabetic retinopathy, retinovascular occlusive disease, and choroidal neovascularization as occurs most commonly in age-related macular degeneration. While highly valuable, fluorescein angiography is invasive, time consuming, and sometimes associated with undesirable side effects [1,2]. Development of a non-invasive, less costly alternative to fluorescein angiography is desired.
Optical coherence tomography (OCT) [3,4] is a technique that has demonstrated its ability to non-invasively image cross sectional structure of the retina. Its improved variation called spectral or Fourier domain optical coherence tomography (FD-OCT) [5,6], due to improved sensitivity and acquisition time [7], allowed for volumetric retinal imaging in clinical setting [8,9]. Functionality improvements such as Doppler OCT and similar methods enable OCT to visualize regions of flow in the retina, but visualization has typically been limited to the major veins and arteries within the eye [10][11][12][13][14]. We have developed a novel motion contrast technique for FD-OCT capable of visualizing a wider range of velocities and flow orientations than Doppler OCT techniques. This method has been demonstrated for in vivo by visualization of vasculature in the zebrafish and the retinal and choroidal vasculature of the mouse [15,16]. In this paper we present the initial results of adapting this motion contrast technique to normal human retinal OCT imaging, demonstrating the microvascular imaging capabilities for multiple layers within the retina.

Theory
The phase information available within Fourier domain optical coherence tomography methods allows for the computation of motion contrast. The most common form of this contrast is Doppler OCT, which calculates the axial flow component of scatterers from the average phase change over time: For this calculation, λ is the central wavelength of the light source, n is the average refractive index of the sample and T is the time separation between phase measurements. Typically, the time separation between phase measurements is equal to the A-scan acquisition time of the system, τ. The maximum resolvable flow with this method is typically limited by the 2π ambiguity of phase measurements called phase wrapping, limiting phase changes between π and -π. For a system with a center wavelength of λ = 855 nm, imaging a tissue with refractive index 1.38 with a time separation between phase measurements of 40 µs, the maximum resolvable flow rate without applying additional processing such as phase unwrapping would be v axial,max (z) = ± λ/4nτ = ± 4.0 mm/s. The minimum resolvable flow with this method is determined by the SNR-limited noise σ SNR (z) for each of the phase measurements [17][18][19]. By limiting phase contrast to OCT signals of greater than 10 dB, the velocity noise level for the image is estimated as v axial,min (z) = ± σ SNR (z)/4nτ = ± 0.4 mm/s. The Doppler OCT approach to visualizing vasculature is limited in several aspects. In practice, Doppler OCT rarely achieves the SNR-limited noise limit due to factors such as transverse scanning while acquiring phase changes over inhomogeneities and texture pattern artifacts in the sample [17,20]. This limits the visualization of velocities which have slow flows or are oriented transverse to the imaging direction. As OCT A-scan acquisition speeds continue to increase, this lower velocity limit will continue to increase and visualization of slow flows will further decrease with this method.
A phase sensitive approach using the variance of phase changes has the potential to visualize a larger dynamic range of flow velocities than Doppler OCT, independent of orientation. Variance motion contrast has been previously developed for OCT and FD-OCT systems, using sequential pixel or A-scan measurements to calculate contrast [21,14,22]. Our motion contrast method, which has previously been demonstrated for microscopy and mouse retinas, utilizes phase measurements between sequential B-scans to achieve the imaging of an entire frame and a much larger time separation [15,16]. By increasing the time separation between phase measurements such as with sequential B-scans, the sensitivity to slower flows is increased significantly while faster flows maintain high variance contrast which can saturate due to phase wrapping. While quantitative flow information may be lost for scatterers that experience significant phase wrapping for this time separation, the variance contrast ensures the visualization of these regions.
The motion contrast method in this paper differs from previous approaches through the implementation of a small transverse displacement between sequential B-scans, resulting in a small spatial separation as well as the time separation between phase measurements. Although the transverse separation between phase measurements will create additional phase noise due to retinal inhomogeneities, the effect will remain comparable to SNR-limited phase noise if the motion is only a small fraction of the transverse resolution [17]. The phase noise estimate for a transverse motion of 10% of the transverse resolution is approximately equivalent to the SNR phase noise from a 10 dB signal. A sliding window of 10 sequential B-scans was chosen for calculating intensity averages and phase variance from the data sets presented in this paper. The statistics chosen are a balance between the noise fluctuations within the contrast image and the transverse imaging resolution.
For an example time separation of 10 ms, the minimum calculated variance of the axial velocity required to match the SNR-limited pure noise case of SNR = 1 is: σ v,axial,min (z) = σ SNR (z)λ/4πnT = 5 µm/s. The variance can arise from multiple sources such as variations in the measured flow component, phase noise caused by transverse sample motion relative to the imaging direction [17], or by non-flow motion such as Brownian motion [15].
Increasing the time separation between phase measurements also increases the phase noise contributions due to bulk sample motion. Bulk transverse motion and errors in axial motion compensation account for the dominant sources of phase noise, depending on imaging parameters chosen and the properties of the sample motion.

Methods
The FD-OCT system used to acquire data sets presented in this manuscript was built at UC Davis. The current configuration uses a spectrometer line exposure time of 40 µs, resulting in an A-scan line rate of 25 kHz. The main components of this instrument have been already described in previous reports [23]. The source is a superluminescent diode (855@75 nm) from Superlum (Carrigtwohill, Ireland) allowing 4.5 µm axial resolution at the retina (n=1.38). The lateral resolution of the system at the retina is estimated at 10-15 µm depending on subject's ocular aberrations. Measured power incident on the subject's cornea was equal to 700 µW. System sensitivity is estimated at ~100dB. With our current spectrometer design, the maximum axial range (seen after the Fourier transform) is 2.7 mm in free space and around 2 mm in the eye.
Data were acquired from three subjects ranging in age from 45 to 57 years. All had normal ocular media and were free of any disease of the retina and optic nerve upon clinical examination at the UC Davis Department of Ophthalmology & Vision Science. Each 3D data set presented in this paper is composed of 1000 B-scans, each consisting of 125 A-scans acquired over 1 mm scanned uni-directionally with an acquisition time of 5.0 ms and a flyback time of approximately 1.6 ms, resulting in a total acquisition time of approximately 6.6 seconds. The transverse spacing between successive B-scans varied for different motion contrast acquisition protocols between 0.5 µm and 1.0 µm. The choice of different sampling patterns for a given subject only altered the spacing between sampling spots and did not affect the transverse resolution. According to the sampling theory, a minimum spacing of 5 µm between sampling points is required to achieve system's nominal lateral resolution. Data sets were selected from the subjects for this paper based on the lack of significant eye motion presenting during the acquisition time.
The flow chart in Fig. 1 describes the data processing operations implemented to produce the motion contrast images from the series of successive B-scans. For the images presented in this paper, alignment between B-scans was performed only in the axial direction z, using the cross-correlation of the intensity images to a selected reference frame. Due to the short acquisition time of each B-scan, motion was assumed to be small enough that no A-scan realignment was performed within each B-scan. The details on the bulk motion removal have been previously described using phase changes between successive aligned B-scans [15]. Following bulk motion removal, no phase unwrapping algorithm was applied to the phase change data. A sliding window of 10 sequential aligned B-scans was used to calculate average intensity images and phase variance contrast images from the data set. Thresholding was applied to the phase contrast data to limit the visualization of scatterers with an average intensity less than 10 dB above the noise level of the image. Additional processing of the three-dimensional contrast data sets primarily consisted of smoothing with a 3D Gaussian filter of 4 µm radius in the transverse directions and 2 µm radius in the axial direction.
Visualization of the three dimensional data sets is presented through summation images, two dimensional images created through a selected depth region of the data sets. Summation images are presented in either fixed thickness summation slices or summation over variable layer thicknesses based on the thickness of specific tissue layers. Edge detection of the RPE layer from each averaged intensity B-scan was used to separate the retinal and choroidal regions of the image, allowing for layer-specific summations to be produced. Additional noise removal has been implemented to some of the data sets to compensate for phase noise created by transverse eye motion during the acquisition time. Since the acquisition time for each B-scan is small, it is assumed that all motion occurring during Bscans is linear and contributed noise uniformly to the phase change measurements of the entire B-scan. With the retinal vasculature expected to comprise only a small fraction of the retinal area anterior to the RPE, this region was used to estimate and remove the phase noise contribution. The values of the phase variance image which were anterior to the RPE and had average intensities greater than 10dB above the noise level were collected in a histogram to estimate the phase noise. The maximum count of each phase variance histogram was used as the bulk phase error estimate and removed from the corresponding phase variance image. The number of bins used for the histograms was manually adjusted for each data set to reduce motion artifacts. The exception to this noise removal method was for B-scans which scanned along the length of a vascular event, which causes the majority of the analyzed region to contain vasculature or artifacts below the vascular locations which affects the noise estimate. To account for this effect, thresholds were applied to remove high variance contrast data from the estimate analysis for cases where a large percentage of the image area contained high contrast associated with large vasculature flow and motion artifacts below those flow areas. Due to the negligible change observed with the removal of expected SNR-limited phase noise for the data processing thresholds used, the phase noise contribution from this source has not been removed in the images presented.

Results and discussion
Cross-sectional image examples of the processed OCT intensity and phase information for several selected stages of the motion contrast data processing are presented in Fig. 2. OCT intensity and phase information available from a single B-scan acquisition are shown in Figs. 2 (a) and (b), demonstrating the retinal structure within the intensity image and a lack of any retinal information within this form of the phase image. A set of 10 sequential intensity Bscans were aligned and the linear intensity was averaged to produce the image in Fig. 2 (c). The decision to perform axial alignment only on the entire B-scans as opposed to the individual A-scan realignment did not produce any significant artifacts within the averaged images. For increased B-scan acquisition times from larger retinal scans, bulk motion variations may require individual A-scan axial realignment and transverse B-scan alignment to maintain the quality of the averaged intensity image.
The effects of bulk motion are observed in the Fig. 2 (d), an image calculated from the phase change between two successive B-scans. The phase-wrapped effects of bulk motion are apparent on the retinal areas of the image; regions without OCT intensity information from the retina appear as a random distribution of phase changes. The bulk motion is observed as vertical lines across the image which changes in value across the image based upon the variations in axial motion during acquisition. Bulk axial motion removal based on a previously detailed method significantly reduces the artifacts within the image, as shown in the example of Fig. 2 (e). Through the phase variance calculation of phase changes across 10 sequential B-scans corrected for bulk motion and the application of a threshold based on the averaged intensity image, the vascular features of the retina become resolvable. Three dimensional smoothing of the data can improve this visualization by reducing specular phase noise appearing within the image while maintaining the smallest motion contrast regions observed.  The data presented in Fig. 3 demonstrated the coarse depth separation of motion contrast between the retina and the choroidal regions, but finer depth separation can easily be achieved. With the depth resolution available from OCT, we are able to segment the vasculature from different depths with a resolution limited primarily by the vascular feature sizes. Within the motion contrast data set used to create the summation images in Fig. 3, we have observed three separate layers of microvasculature, located at depths of approximately 60 µm, 100 µm, and 150 µm below the anterior surface of the retina. These correspond to depth regions associated with the posterior part of the ganglion cell layer (GCL), the anterior region of the inner nuclear layer (INL), and the posterior region of the INL within the retina. Segmentation of these individual layers is presented in Fig. 4.
To optimize the visualization of each vasculature layer, a summation image was created for each layer summing over a 16 µm thick depth region centered at each depth and presented in Fig. 4. This is presented alongside example B-scans of the intensity and motion contrast used to register the depths used for these three separate image layers. As demonstrated in previous work [16], artifact noise of motion contrast appears below strong flow regions and is only reduced or removed for depth regions where the intensity signal is reduced below threshold levels. For this reason, the first two anterior layers appear to have the same major vasculature and some microvasculature differences. To aid in visualizing the different retinal vascular layers, a composite image of the vasculature observed in the three retinal layers is presented in Fig. 4. In the motion contrast image, the average of the two vascular layers in the GCL and the anterior of the INL are presented in green and overlaid with the posterior INL layer in red. As demonstrated in Fig. 4(f), the top two vascular layers observed shows minimal spatial correlation to the vascularization in the posterior INL, confirming that motion contrast regions in this image are not artifacts from microvascularization anterior in those retinal locations. The depth sensitivity of the microvascular imaging is a feature that cannot be replicated by clinical fluorescein angiography. Reducing the transverse sampling by a factor of two in the data set acquisition does not degrade the vascular visualization capabilities with this technique, as demonstrated in Fig. 5. A retinal scan area of 1 mm x 1 mm was acquired in the same acquisition time as the previously shown data set at a similar retinal location. As expected, only the major retinal vasculature is visible within the intensity summation image based upon the blood absorption. Within the retinal depth summation image of the motion contrast data, the major vasculature is observed as well as the network of motion contrast which appears to be microvasculature. To our knowledge, the level of detail observed within this retinal vasculature has not been duplicated within other phase sensitive OCT imaging techniques. Additional noise removal from the motion contrast data sets is critical to maintain vascular visualization over the entire retinal scan area. While microsaccadic motion can be easily identified within the acquisitions, transverse eye motion caused by slower changes in fixation position also contributes to the noise observed and can be much more difficult to compensate. Fluctuations in this type of transverse eye motion during acquisition can cause a varying level of phase contrast noise across the entire scan area. Due to the relatively short acquisition time of an individual B-scan relative to the expected motion drift variations, the noise experienced due to motion drift is approximated as a constant contribution across each phase variance Bscan image. Histogram-based analysis of the motion contrast over the retinal portion of the image estimated the noise contribution, assuming that the majority of the cross-sectional image contains static retinal layers as opposed to retinal vascular regions and motion artifacts. The motion contrast summation images demonstrated in this paper have implemented this noise removal technique to maximize vascular visualization.
Applying this motion contrast imaging technique to the foveal region of the retina demonstrates the extent of the vascular imaging capabilities, since the fovea contains a significant avascular region. Within this region, there is expected to be no motion contrast within the retinal depth. Figure 6 (b) -(d) contains 1 mm x 1 mm retinal data acquired over the foveal region for motion contrast calculations. The 3 mm x 3 mm intensity summation image presented in Fig. 6 (a) was used only for context to identify the approximate foveal region imaged and does not have any associated motion contrast information. The motion contrast observed in Fig. 6(d) provides a valuable tool to visualize the avascular region boundary within the fovea. While the major retinal vasculature which can be observed within the intensity images only extends over the right half of the imaged region, the motion contrast images demonstrate a network of microvasculature extending further towards the center of the fovea. The edge of the observed microvasculature structure is associated with the initial thinning of the retina within the OCT images at the fovea. Figure 6(c) demonstrates the retinal motion contrast summation image before the additional noise removal. Horizontal black lines are associated with contrast data which was eliminated from the image due to thresholds imposed on the processing of the motion contrast data. Each of these lines is associated with a microsaccade occurring during the image acquisition. The horizontal lines on the image which appear to have high contrast have been identified as motion contrast artifacts which have not been compensated with bulk axial motion noise removal and thresholding. The motion contrast image after additional noise processing removes almost all of these artifacts and improves vascular visualization across the entire image. Improvements to the demonstrated motion contrast technique focus around increasing the retinal scan area while maintaining or improving vascular visualization capabilities. Due to the fixation limitations of typical subjects, acquisition times should remain within a few seconds or retinal tracking is required. This contrast method can be adapted to faster OCT systems without significant change in vascular visualization, allowing for larger retinal scan areas in the same acquisition time. Further improvements to increasing the scan area can be achieved through stitching together multiple data acquisitions guided by retinal tracking information or retinal features [26]. Figure 7 presents the comparison of the motion contrast images with a fluorescein angiography image of the same subject cropped over a selected region of interest. A retinal area of the fluorescein angiography image of approximately 2 mm x 2 mm was chosen which includes the OCT data scan regions corresponding to the data presented in Figs. 3 and 5. Figures 3 (b) and 5 (c) were overlaid over the fluorescein angiography image based on the approximate scan location. Distortions caused by transverse motion during the acquisitions were not compensated for in this comparison, which accounts for position discrepancies of the major vascular features between the images. Comparison of microvascular imaging between the OCT motion contrast images and the fluorescein angiography image shows comparable vascular imaging capabilities, with similar feature sizes and spacing between vasculature within the images. Direct comparison between specific microvascular events is difficult due to distortions in the OCT image which limit the accuracy of the vascular location, and by the contrast limitations within regions of the angiography image but the visualization capabilities appear to be analogous. Detailed analysis in comparing the microvascular visualization capabilities between fluorescein angiography and the demonstrated phase variance motion contrast technique requires an accurate mapping of retinal locations, currently limited by transverse motion during data acquisition. Retinal tracking and faster data acquisitions are both options available to reduce the effect of these motions and improve mapping of the vasculature. The demonstrated motion contrast images in this paper have shown resolution-limited features independent of vascular orientation or flow velocity, and with improvements to the retinal scan area while maintaining this visualization, this technique is expected to continue to produce comparable vascular visualizations to fluorescein angiography over larger retinal scan areas.

Conclusion
Phase variance-based motion contrast has been demonstrated for several different regions of in vivo human retina. Microvasculature visualization in the retina with this technique was shown to have a capability superior to other techniques previously demonstrated with OCT imaging. This motion contrast technique has been demonstrated to be analogous to fluorescein angiography with the ability to visualize resolution-limited vasculature independent of vessel orientation and flow velocity. Depth segmentation of the motion contrast data has demonstrated three separate vascular layers within the retina, a feature unparalleled with clinical fluorescein angiography. Foveal imaging has demonstrated microvasculature visualization extending further into the avascular region beyond the end of the major retinal vasculature. Further improvements to the motion contrast imaging technique will improve vascular visualization, reduce the effects of noise attributed to fixation motion, and increase the retinal scan area of the images.