Comparison of amplitude-decorrelation, speckle-variance and phase-variance OCT angiography methods for imaging the human retina and choroid

: We compared the performance of three OCT angiography (OCTA) methods: speckle variance, amplitude decorrelation and phase variance for imaging of the human retina and choroid. Two averaging methods, split spectrum and volume averaging, were compared to assess the quality of the OCTA vascular images. All data were acquired using a swept-source OCT system at 1040 nm central wavelength, operating at 100,000 A-scans/s. We performed a quantitative comparison using a contrast-to-noise (CNR) metric to assess the capability of the three methods to visualize the choriocapillaris layer. For evaluation of the static tissue noise suppression in OCTA images we proposed to calculate CNR between the photoreceptor/RPE complex and the choriocapillaris layer. Finally, we demonstrated that implementation of intensity-based OCT imaging and OCT angiography methods allows for visualization of retinal and choroidal vascular layers known from anatomic studies in retinal preparations. OCT projection imaging of data flattened to selected retinal layers was implemented to visualize retinal and choroidal vasculature. User guided vessel tracing was applied to segment the retinal vasculature. The results were visualized in a form of a skeletonized 3D model.


Introduction
In human ocular anatomy two distinct vascular systems are recognized in the histology and angiographic imaging of the fundus [1,2]: the retinal vasculature and the choroid. Both systems are derived from the ophthalmic artery. The retinal vasculature is responsible for maintaining the proper function of the inner retinal layers (from the nerve fiber layer to the outer plexiform layer) while the choroid sustains the function of the photoreceptor and retinal pigmented epithelium (RPE), and regulates the thermal conditions of the ocular tissues. A range of eye diseases can affect the vascular layers. Most notably, the two most prevalent diseases: diabetic retinopathy and age-related macular degeneration (AMD) compromise the retinal vasculature and the choroid, respectively. Therefore, the capability of non-invasive imaging of retinal and choroidal vasculature and assessment of their integrity is an important asset in clinical diagnostics. The retinal vasculature emerges from the retinal artery located in the optic nerve head and covers the retina with a network of arteries and veins located in the superficial layers, close to the nerve fiber layer. The larger retinal vessels divide into a network of smaller vessels down to the capillary level. The first and most anterior capillary system is located in the nerve fiber layer. The second capillary system is formed by vessels ramifying towards the outer retinal layers and is located in the inner plexiform layer (the inner capillary plexus). The third capillary system is located in the outer plexiform layer (the outer capillary plexus). The nuclear layers are largely devoid of vessels, with the exception of the inner nuclear layer which contains vessels supplying and draining the outer capillary plexus. In the healthy retina there are no vessels in the outer retinal layers: the outer nuclear layer, photoreceptor layers and in the retinal pigment epithelium (RPE). The metabolites necessary to sustain the function of these structures are provided by the choroidal vasculature by way of diffusion through Bruch's membrane from below.
The choroid is located between Bruch's membrane and the sclera. The choroidal vasculature originates from the ciliary arteries penetrating through the sclera. The network of large choroidal vessels (Haller's layer) divides into a system of medium size vessels (Sattler's layer). Sattler's layer supplies the thin layer of choriocapllaris located underneath Bruch's membrane. Choriocapillaris is a meshwork of densely packed, interconnected vessels forming characteristic mesh-like structures. It is organized in a lobular structure [3]. Blood with nutrients is delivered by arterioles located in the centers of the lobules. Blood with metabolic waste materials is collected by peripheral venules surrounding the lobules.
Both vascular structures of the eye fundus, the retinal vessels and the choroid, are critical in maintaining the proper function of the retina. Methods allowing for their visualization are important in monitoring the health of ocular tissues and in disease diagnostics. The most commonly used vascular diagnostic methods in ophthalmology are fluorescein angiography (FA) and indocyanine green (ICG) angiography. In recent years advances in optical coherence tomography (OCT) imaging have led to development of optical coherence angiography (OCTA) methods which are providing additional, clinically relevant information about vascular structures of the eye.
OCT angiography is a group of optical coherence tomography imaging methods facilitating visualization of blood flow in biological tissues. Similar to Doppler OCT techniques, OCTA detects motion in the imaged sample. The OCT angiography methods, however, ignore the quantitative measurement of blood flow and use motion as a contrast mechanism to visualize the location of moving cells. OCTA relies on interferometric detection of low temporal coherence light back-scattered by moving blood cells. The interferometric detection is usually realized in a Michelson or a Mach-Zehnder interferometer configuration. In Fourier-domain OCT methods, the signal is derived from a spectral fringe pattern recorded on either a spectrometer (as in spectral domain techniques) or with a photoreceiver during wavelength sweeps of specially designed light sources (swept-source techniques). Its Fourier transformation yields the OCT signal in the form of complex amplitude (phasor): which describes the amplitude A 0 (x,y,z;t) and phase Φ(x,y,z;t) of the spectral interference fringe pattern corresponding to each location (x,y,z) of the imaged object, at each time point t.
The OCT measurement provides discrete sampling of A(x,y,z;t): v = (p, q, l) -indexes of the OCT data voxels, t m -discrete time points. When the blood cells flow across the imaging beam, they cause fluctuations in the registered interference signal (raw data) which, after Fourier transformation, give rise to depth-localized fluctuations in the resultant complex amplitudes. To detect motion, a series of OCT data has to be acquired at the same location of the object at time intervals Δt = t m + 1 -t m sufficient to observe changes of the complex amplitude introduced by the flow. The design of scan protocols must take into account flow velocities expected in the imaged object to provide correct timing of the data acquisition for the OCTA analysis [4][5][6]. In the currently available research grade OCT systems operating at A-scan acquisition rates up to 100kHz (100 000 Ascans/s), the most commonly used OCTA data acquisition regime is MB-scan mode, in which a series of M subsequent B-scans is acquired at the same location of the sample. The time interval Δt between data points taken for OCTA analysis at a given voxel v is determined by the acquisition time of a single B-scan. Thus, Δt is given by the B-scan numbers m in the MBscan, and the complex amplitude at each voxel v can be described as: To detect and visualize fluctuations of the OCT signal, various metrics can be defined giving rise to different OCT angiography methods. A wide range of OCTA techniques has been developed by different research groups resulting in a vast literature on this topic . This provides ample choice of methods for application in clinical imaging. But it also raises questions, and often confusion, as to the differences between results obtained with various OCTA modalities and which of them would be the best choice in imaging of specific vascular systems of the eye. Only a few papers have included comparisons of different methods [34, 44, 46, 47] and even fewer papers report quantitative comparisons of OCTA performance in the imaging of the choroid.
In this paper we focus on empirical comparison of three OCTA methods: a) speckle variance (SV OCTA), b) amplitude decorrelation (AD OCTA), c) phase variance (PV OCTA). This collection is not exhaustive but does includes the techniques most commonly discussed and implemented for ophthalmic applications.
The first two methods analyze amplitude changes of the OCT signal to visualize the flow and ignore the information contained in the phase. The third relies on the detection of phase changes but ignores the amplitude variations.
In the speckle variance method, the variance of the amplitude fluctuations between Bscans is calculated to visualize flow [17,22] as defined by: where M is the number of B-scans in the acquired MB-scan series. In the locations where the flow is present, the variance is higher than in the signal originating from static tissue. The amplitude-decorrelation method uses the correlation as a metric to detect the changes in the OCT signal. The correlation function of the discrete, complex, OCT signal (Eq. (2)) can be defined as: Although the complex correlation can be used to visualize flow [24], in the present study we follow the approach in which only the amplitude of the OCT signal is used to calculate the decorrelation [34]: which can be also written as (Eq. (3)): In the areas of flow, the decorrelation of the amplitude calculated between the B-scans is higher than in the static tissue, where it is caused by the noise. A high decorrelation value, therefore, implies the presence of vessels in the imaged tissue.
The phase-variance method uses a variance of the phase differences Δφ v,m occurring between B-scans as a metric to visualize flow [8,12,20]: , The phase differences used in the phase-variance method correspond to the phases of the complex correlation function (Eq. (5)). In the areas of flow the variance is higher than in the areas of static tissue and indicates the presence of vessels. All OCTA methods can benefit from application of averaging methods to improve the signal-to-noise ratio and in consequence to improve the visualization of the vasculature. In this paper we explore two averaging methods: a) split spectrum, b) volume averaging. The split-spectrum method was developed in connection with the amplitude decorrelation OCTA method giving rise to split-spectrum amplitude-decorrelation angiography (SSADA) [34,48]. It relies on splitting the acquired interference spectra into narrower, usually overlapping bands. The splitting is usually performed by windowing the spectra with a Gaussian function. Each separate band is Fourier transformed to obtain a set of A-scans. Decorrelation B-scans are then calculated for each band and averaged to generate the final OCTA images. Although the split-spectrum method results in images of lower axial resolution than the original OCT data, it was shown to improve the visualization of vasculature in the human retina [34,46].
Spatial averaging methods were introduced to the OCT imaging techniques to improve the signal-to-noise performance in the tomograms, to reduce speckle noise and to reduce eye motion artifacts [49-57]. Volume averaging relies on acquisition of several 3D data sets and averaging of the corresponding OCTA volumes. In current OCT systems, the volumetric data sets are usually acquired in 3-6 s. During this time involuntary eye motions may introduce artifacts in the acquired data. Careful design of scan protocols and data processing methods is, therefore, recommended to reduce the artifacts and to avoid the loss of spatial imaging resolution. Several volume averaging methods with motion correction have been introduced to OCT and OCTA imaging by different research groups. Some groups have performed post processed motion correction based on acquired OCT data [56-58] while others rely on active (real time) eye tracking [59][60][61][62][63]. In this study we have implemented motion correction and volume averaging in data post processing, using a method similar to [57].
We compared the performance of three OCTA techniques (speckle variance, amplitude decorrelation and phase variance) for imaging of the human retina and choroid. Two averaging methods were used (split spectrum and volume averaging) to improve visualization of the vascular structures. The OCTA techniques along with averaging methods are compared quantitatively in their visualization of the choriocapillaris layer. Although several metrics have been proposed for evaluation of vessel visualization in OCTA images [22,34,46, 64], they were primarily designed for retinal vasculature imaging and are calculated from en face OCTA projections. These approaches are not well suited to evaluate the visualization of the choroid. We proposed to calculate contrast-to-noise ratio (CNR) between the photoreceptor/RPE complex and the choriocapillaris layer as a metric for evaluation of the OCTA performance in static tissue noise reduction and imaging of the choroid. We demonstrated that implementation of intensity based OCT imaging and OCT angiography methods allows for visualization of retinal and choroidal vascular layers known from anatomic studies with retinal preparations [1]. OCT projection imaging of data flattened to selected retinal layers was implemented to showcase retinal and choroidal vasculature Userguided vessel tracing was applied to segment the retinal vasculature. The results were visualized in the form of a skeletonized 3D model.

Study participants
Imaging of human subjects was performed under a protocol approved by the UC Davis Institutional Review Board. Three normal subjects (N1, N2, N3, ages: 63, 35, and 29 years), a patient diagnosed with geographic atrophy (GA, age 66), and a patient with age-related macular degeneration (AMD, age 64) were imaged with a swept-source OCT system. Subject preparation included instillation of eye drops: 1% Tropicamide and 2.5% phenylephrine for pupil dilation and cycloplegia. During imaging, head position was stabilized with a forehead rest and bite bar. The imaged retinal location was selected by guiding the gaze with a white fixation cross displayed on a liquid crystal display (LCD). Light exposures were at the level of 1.3mW, which is below the maximum of the ANSI laser safety standards [65, 66]. OCTA volumes were acquired with a swept-source OCT system (Fig. 1). The light at 1040nm central wavelength was generated by a swept-source laser (Axsun Technologies, USA) operating at 100 000 sweeps/s. An 80:20 fiber coupler was used in the Mach-Zehnder interferometer to split the light between the object arm (20% channel) and the reference arm (80% channel). A standard OCT design was used for the construction of the object arm, which consists of: a collimator L c (f c = 11mm), XY galvanometer scanner, and two lenses (L 1 , L 2 , effective focal lengths: 100 mm and 75 mm). The light returning from the object arm was directed to the detection system via the 80% transmission channel of the fiber coupler. A cube retro-reflector was used in the reference arm to direct the light from the 80:20 coupler via an open air path to the collecting fiber and further to the detection system. An iris with controllable diameter was used to adjust the light power in the reference arm. The reference arm fiber was inserted into a polarization controller to match the polarization of light propagating in the object arm. Dispersion compensation was performed by incorporating a fiber of an experimentally determined length to the reference arm. The residual dispersion mismatch was corrected in data post processing. The system specifications of the OCT setup are summarized in Table 1.

Swept-source OCT system
The OCT signal was acquired with a dual-balanced detection system containing a 50:50 fiber coupler and a pair of photodiodes (Thorlabs, USA). A data acquisition card (Alazar Technologies Inc., Canada, ATS9350, 12-bit resolution, of 500MS/s sampling rate,) and a graphic processing unit (GPU) based data acquisition software (JXW Software Ltd, Canada [67-69]) were used for data collection. To account for the sweep to sweep jitter of the sweptsource laser [70, 71], a fiber Bragg grating (FBG) was used [72] in one of the 50:50 detection channels. The FBG introduced a notch to the registered spectra at the beginning of each sweep at 989 nm. The notches were used to correct the sweep position relative to the beginning of each acquisition window. Prior to OCT data processing, the parts of the spectra containing the notches were removed.

OCT angiography scan protocols
Multiple B-scans (MB-scans) were acquired at each retinal position of the slow scanner enabling detection of the blood flow induced OCT signal changes. The details of scan protocols are shown in Table 2. In their design we avoided undersampling (i.e., the step size was smaller than half the spot size of the beam), kept the imaging time to a minimum, and limited volume sizes to ~1GB. The visualization of Haller's layer in the choroid was achieved with a wide-field scan in intensity only OCT imaging. The details of the scan protocols are presented in Table 3.

OCT data processing
The acquired OCT data were linear in wavenumber space, therefore Fourier transformation could be performed directly with the collected spectra. However, numerical corrections were applied prior to transformation to account for the sweep to sweep jitter of the laser (as described in paragraph 2.2.), to compensate the residual dispersion mismatch between the object and reference arm of the interferometer (dispersion compensation), and to correct the non-Gaussian shape of the laser emission spectrum (spectral shaping). Amplitude and phase were extracted from the Fourier transformed data for further OCT angiography analysis.

OCT angiography data processing of B-scans
In speckle variance OCT (SV OCTA) and amplitude decorrelation (AD OCTA) methods, only the amplitude of the complex signal is utilized to visualize the flow. To generate the SV OCT angiography images we used Eq. (4) and implemented data processing similar to [22].
To realize the AD OCTA we used Eq. (7) and based our data processing on [34]. In our implementation, the processing steps used for both methods (SV and AD OCTA) are as follows (details are given in Appendix 1): a. Thresholding of the intensity images to remove intensity noise which would introduce erroneously high variance or high decorrelation values.
b. Correlation of the B-scans within the MB-scan series to correct for the shifts of the intensity distribution (speckle pattern) caused by eye motion and mechanical instabilities of the OCT system.
c. Calculation of the variance or decorrelation of the amplitude, depending on the OCTA method.
d. Removal of the noisiest component.
In the phase variance OCT method (PV OCTA) only the phase of the complex signal is used to compute the flow signal. Our implementation of the PV-OCTA method was based on the data processing described in [12]. The following data processing steps were implemented.
a. Intensity image thresholding to remove the noise pixels before the phase data processing.
b. Bulk motion removal based on the averaged shifted histogram analysis of the phase data [10].
c. Calculation of the phase variance according to Eq. (8).
In each of the OCTA data processing methods, average intensity B-scans were also calculated from the sets of B-scans in the MB-scan series. The average B-scans were used in further image processing and data visualization procedures.

Implementation of the split-spectrum method
The split-spectrum method was originally introduced as an enhancement of the amplitudedecorrelation method. In our study we extended its application to the speckle-variance and phase-variance methods. In the spectrum splitting procedure we followed two arbitrary rules of thumb. We kept the axial imaging resolution at the level approximately matching the sizes of the smallest imaged features (capillary vessels). And, following [34], we kept the axial and transverse imaging resolutions the same. The spectra were split into four equally spaced, partially overlapping bands by Gaussian windowing of the acquired interferometric signals. The full width at half maximum (FWHM) of the Gaussian window was set to be half of the FWHM of the emission spectrum of the light source decreasing the axial imaging resolution to ~12 µm. OCT angiography computations were performed for each split spectrum and the resultant data averaged to yield OCTA images. Further data processing steps were identical to those implemented with the full interference spectra, and are described in the following sections.

Image processing and visualization
Implementation of the OCT angiography data processing generated sets of OCTA cross sections and corresponding sets of average intensity B-scans. These data were visualized in the form of OCT en face projections by axial summation with Gaussian weighting of depth slices selected from 3-D data sets. Images obtained with all three OCTA methods underwent identical image processing steps to provide a common baseline for further comparisons. Image processing and visualization consisted of the following steps (details are given in Appendix 2): a. Removal of bulk motion affected cross-sections by automatic detection of OCTA Bscans with the mean signal intensity exceeding an empirically selected threshold. The locations of motion affected cross-sections were stored as a separate vector to facilitate transverse motion correction and volume averaging procedure described in section 2.8.
b. Axial registration of intensity B-scans was performed to correct for axial (depth) displacements of the cross-sections caused by the eye motion. All B-scans were cross correlated with the center B-scan in the volume data set. Only the axial component of the displacement was corrected. Transverse displacements were ignored at this data processing step.
c. Flattening of the retinal images to selected anatomical layers was performed to allow for visualization of different vascular layers in the retina and choroid. Flattening to the inner limiting membrane was performed to visualize the retinal vasculature. Visualization of the outer capillary plexus was performed in the data flattened to the outer plexiform layer. The projection images of the choroid were generated from the data flattened to the RPE.

Motion correction and volume averaging of the OCT data
The transverse motion correction and volume averaging methods used in this study are similar to a concept published in [57]. Our method was developed under simplifying assumptions that: -Data blocks between the saccadic motions of the eye are free of motion artifacts, -Saccadic motion causes only translations of the data blocks. Identical motion correction and volume averaging was applied to images obtained with all compared OCTA methods. Four OCT data sets were acquired at the same location in the retina. A user-supervised registration procedure was applied to the data flattened to the RPE. The transverse motion correction was performed in projection images of the inner retinal vasculature. Each volume was broken into sub-volumes at the locations of the removed motion affected B-scans. En face projection images of overlapping sub-volumes were overlaid with each other by cross-correlation, providing transverse shifts for motion correction. The sub-volumes with corrected transverse positions were registered axially and averaged. Further details of the method are given in Appendix 2.

Average OCT angiography depth profiles and normalization of the OCTA data
A mean of all OCT angiography A-scans ( Fig. 2(g)) was computed from the volumes flattened to the RPE. The average depth profiles facilitated comparison of the OCTA methods, and provided guidance in selection of depth locations for projection image generation, as well as for OCTA data thresholding and normalization.
In the data thresholding and normalization procedure we took advantage of the anatomical structure of the retina and choroid. The outer nuclear layer (ONL) is a highly transparent tissue devoid of vasculature and therefore should not produce OCT angiography signal. The residual signal registered in this layer was considered as noise and used to threshold the OCTA data. The mean and standard deviation of the signal from the ONL were calculated. Voxels with signal values smaller than three standard deviations above the outer nuclear layer mean were set to zero. In contrast, the OCT angiography signal from the choriocapillaris produced the strongest peak in the average depth profiles due to the blood flow in a dense vascular meshwork. The maximum value of the OCTA signal in the choriocapillaris layer was used to normalize the OCTA data sets. The thresholding and normalization procedures, identical for all three of the OCT angiography techniques, were implemented to provide a common baseline to compare the techniques.

Generation of en face projections
Gaussian windowing was applied to the A-scans to select the layers of interest and generate en face projections of the retinal and choroidal layers. A Gaussian function provided weighting coefficients in the axial summation of intensity to gradually decrease the influence of pixels located away from the center of the window on the resultant image. The windows were constructed in such a way that the full width at half maximum (FWHM) of the Gaussian function corresponded to 1/3 of the user selected depth range. All points outside the user selected range were clipped to zero. En face projections were generated by depth summation of the windowed data volumes. The FWHM thicknesses and locations of the Gaussian windows depended on the visualized vascular layer. Details of the windowing parameters applied to the data obtained from the three normal subjects are given in Table 4. The illustration of Gaussian window positions in the average depth profiles is provided throughout the Results section in the relevant figures.

Contrast-to-noise ratio of the choriocapillaris imaging
The OCT angiography methods should reject OCT signal originating from static tissue while preserving the signal from the moving scatterers (flowing blood cells). The photoreceptor-RPE complex consists of highly scattering tissues (RPE, photoreceptors outer segments tips and photoreceptors inner/outer segments junction -IS/OS) producing strong signal in the intensity images [73]. However, in a healthy eye there are no blood vessels in these layers and no flow generated OCTA signal should be detected. The residual OCTA signal registered in the RPE complex indicates how well static tissue signal is rejected by different angiography methods.
To compare the ability of the OCT angiography methods to reject the static tissue noise without reducing the signal from the flowing blood, the contrast-to-noise ratio (CNR) between the choriocapillaris and the photoreceptor layer including IS/OS junction was calculated as where CC I , PR I are the mean values of the OCTA signal in the choriocapillairs and photoreceptor layers, respectively, and PR σ is the standard deviation of the OCTA signal in the photoreceptor IS/OS layer. The CNR was calculated in the thresholded and normalized OCT angiography volumes.

Overview of the vascular layers
Average OCT angiography profiles aided the selection of depth ranges used for generation of the en face projections by Gaussian depth windowing. Figure 2 shows example profiles obtained in the SV, SSADA and PV OCTA data sets flattened to the Inner Limiting Membrane (ILM), Outer Plexiform Layer (OPL) and the RPE. Peaks corresponding to three layers: inner plexiform layer (IPL), OPL, and choriocapillaris were selected as depth references to visualize different retinal and choroidal vascular systems shown in Figs. 3-7. All major retinal layers containing vessels are visualized as peaks and layers devoid of vessels as valleys. However, highly scattering tissues produce a non-zero signal even if they do not contain vessels. This may be considered OCTA noise. An example is the photoreceptor-RPE complex visible as a high "plateau" in the average profiles of the data flattened to the RPE. There are the two main contributions to this signal. The first is the signal intensity-dependent OCTA noise [12]. The second is the signal generated by blood flow in the overlying retinal layers (so-called "Doppler traces") [74, 75]. Figures 3-7 provide an overview of vascular layers visualized with the OCTA methods. The images of the retinal vasculature ( Fig. 3-5) were generated from data sets obtained with the SSADA algorithm and application of volume averaging. We selected this method to showcase the capabilities of the OCTA techniques using results of the comparisons shown later in the paper. SSADA with volume averaging provided the clearest images of retinal vasculature. Similarly, split spectrum PV OCTA with volume averaging rendered the most detailed images of the choriocapillaris and Sattler's layer and therefore was selected for the overview of the vasculature in these layers. The intensity en face projections are shown below corresponding OCTA projections to demonstrate the relative visualization improvement of the indicated OCTA technique.
In Fig. 3 major retinal vessels are visualized in the projection which includes three anatomical layers: the ganglion cell layer, the inner plexiform layer and the inner nuclear layer (Fig. 3(a)). The meshwork of the inner capillary plexus is visualized in the inner plexiform layer (Fig. 3(b)). In comparison with corresponding intensity projections (Fig. 3(d) and Fig. 3(e)), the OCTA projections show more detailed images of the vasculature. In Fig. 4 a sharp image of the outer capillary plexus is generated in the OCTA projection of the outer plexiform layer (Fig. 4(a)). However, an image with similar detail can be obtained from the intensity projections by simple intensity thresholding (Fig. 4(c)). Most interestingly, the projection of the inner nuclear layer (Fig. 4(b)) reveals vessels connecting the outer capillary plexus with the inner capillary plexus and major retinal vessels. As an example we show only results obtained with the SSADA method, although these vessels were visible in all discussed methods. In addition, visualization of connecting vessels does not suffer from "Doppler traces/shadows" artifacts because of their location in a highly transparent tissue (INL -inner nuclear layer). This finding suggests that OCTA techniques may allow for retinal vessel tracing in the human eye in vivo down to the level of capillary plexuses. User supervised tracing of the vessels is presented in Fig. 5 Fig. 2(a). Projections of the inner plexiform layer capillaries were obtained by Gaussian depth windowing indicated by the blue line in Fig. 2(a). Projections of the nerve fiber layer capillaries were obtained from full segmentation of the NFL (i.e. axial summation was performed between anterior and posterior boundaries of the NFL). Image location: 6° nasal, 4° inferior relatively to the fovea. Image sizes: 1.8 x 1.8 mm.
The human retina contains yet another layer of capillary vessels located in the nerve fiber layer. The en face projection generated in this layer (Fig. 3(c) and 3(f)) shows fragmented, string-like bright features which may correspond to the nerve fiber layer capillaries. For comparison, the intensity projection shows highly reflective nerve tissue with fragmented outlines of larger vessels, however no evidence of capillaries is visible.
The images of the choroidal vasculature presented in Figs. 6 and 7 were generated with the split-spectrum phase-variance OCTA method along with volume averaging. The most prominent peak in the OCTA profiles in the normal eye (Fig. 6a) indicates the edge of the choroidal complex (choriocapillaris, Sattler's layer and Haller's layer). It has a steep edge toward the RPE side (positive depth values), sharply marking the boundary between the RPE and the choroid. For comparison, an average OCTA depth profile of an eye with geographic atrophy is shown in Fig. 6(d). Although the prominent peak corresponding to the choroid is visible, it lacks the sharp edge on the side of the RPE. An example OCTA cross-section of the normal eye ( Fig. 6(b)) shows a bright band of signal located underneath the RPE, where the choriocapillaris is expected. The cross section obtained in the eye with geographic atrophy (Fig. 6(e)) lacks the bright band under the RPE. Instead, patches of strong OCTA signal are randomly scattered suggesting the presence of larger, more separated vessels. In Fig. 6(h), an OCTA projection located directly below Bruch's membrane of the normal eye shows a bright meshwork of OCTA signal with randomly distributed dark spots suggesting a tight mesh of capillary vessels, with some unwanted intrusion of retinal vessel shadows. The OCTA projection selected at the same location underneath Bruch's membrane in the geographic atrophy case (Fig. 6(j)) shows choroidal vessels (most likely Sattler's and Haller's layers) but   Fig. 2(d).
Projections of the inner nuclear layer vessels were obtained by Gaussian depth windowing indicated by the teal line in Fig. 2(d). Yellow circles in b) indicate selected vessels connecting the outer capillary plexus with the retinal vessels, associated with corresponding clusters of capillaries indicated in a). Image location: 6° nasal, 4° inferior relative to the fovea. Image sizes: 1.8 x 1.8 mm. The results of imaging of choriocapillaris and Sattler's layer in subject N1 are presented in Fig. 7. The choriocapillaris layer has a characteristic mesh-like appearance in the OCTA projection. In comparison, the intensity projection (Fig. 7(c)) lacks the capability of showing any distinct vascular structures, except of the shadows of retinal vessels. In Sattler's layer, the outlines of the vessels are visible in the intensity projection, however the OCTA image reveals more details of their structure. Fig. 6. Comparison of choriocapillaris imaging in a normal eye (top row, subject N3) and a geographic atrophy case (middle row, subject GA) obtained with the split spectrum PV OCT method with volume averaging (Av.). Purple in depth profiles denotes Gaussian depth windowing applied to the data flattened to the RPE. Bottom row: en face projections. N3 image sizes 1.8 mm by 1.3 mm located 8° nasal and 4° inferior from the fovea. GA image sizes: 1.6 mm by 1.6 mm located at 3° nasal and 3° inferior from the fovea.
The visualization of Haller's layer vessels is challenging due to the location of the vessels underneath other highly scattering tissues and due to the high speed of the blood flow. In current OCT imaging systems, Haller's layer vessels do not produce OCT signal and appear as dark patches submerged in the highly scattering tissue. Since the vessels are large, their structure can be better appreciated in wide-field OCT intensity imaging [9,14,33,38,40,41,[78][79][80][81][82][83][84][85][86][87][88][89][90]. In Fig. 8 a comparison of Haller's layer visualization between a normal eye and two pathology cases, AMD and GA, is shown. Although Haller's layer is visible in all eyes, the AMD and the GA case show higher vessel contrast. Most interestingly, in the AMD case ( Fig.  9) a projection located underneath Haller's layer, at the expected level of scleral tissue, reveals circular features which correspond to elongated structures in the B-scans. This suggests the presence of steeply oriented vessels, possibly short ciliary arteries or veins of lymphatics [14,82], which can be better appreciated in Visualization 2 and Visualization 3. Figures 10 and 11 compare retinal vasculature images obtained with the three OCTA methods. The capability of imaging in highly scattering tissue is compared in the inner retinal layers (Fig. 10). More favorable imaging conditions (higher tissue transparency) are found at the edge between the outer plexiform layer and inner nuclear layer where the outer capillary plexus is located. Comparison of the outer capillary plexus imaging is shown in Fig. 11.

Comparison of OCTA methods for imaging of the retinal vasculature
All methods provide similar images of the retinal vasculature. Visual assessment of the OCTA projections suggest that application of the split spectrum method marginally improves the visualization contrast and continuity of vessels. Volume averaging greatly improves the images, particularly in visualizing the outer capillary plexus. High contrast, smooth, continuous vasculature is visible in all three methods. The difference in the imaging quality between the full spectrum and split spectrum is marginal and therefore images from volume averaging of the full spectrum data were omitted.  Fig. 2(g). Projections of Sattler's layer vessels were obtained by Gaussian depth windowing indicated by the red line in Fig. 2(g). Image location: 6° nasal, 4° inferior from the fovea. Image sizes: 1.8 x 1.8 mm.      11. Comparison of performance of OCTA methods in imaging of the outer capillary plexus in subject N1. First column: amplitude-decorrelation method. Second column: phasevariance OCT. Third column: speckle-variance OCT. First row: full spectrum, single volume data. Second row: split spectrum, single volume data. Third row: average (Av.) of 4 volumes obtained with the split spectrum data. Image location: 6° nasal, 4° inferior from the fovea. Image sizes: 1.8 x 1.8 mm.

Comparison of OCTA imaging of the choroid
Comparison of the performance of OCTA methods in choroidal imaging in the normal eye is demonstrated in Figs. 12-15. The first difference between the methods is apparent in the OCTA depth profiles (Figs. 12 and 13). The intensity based methods (AD OCTA and SV OCTA) have lower signal strength originating from the choroid than the phase-based method (PV OCTA). However, application of volume averaging diminishes the difference. The second apparent difference is the height of the "photoreceptor plateau" measured relative to the choriocapillaris peak. The best suppression of the photoreceptor-RPE complex signal (static tissue noise) resulting in the lowest photoreceptor plateau was achieved with the PV OCTA methods, followed by the AD OCTA method and with SV OCTA showing the poorest performance.
To provide a measure of the static tissue signal suppression, the CNR between the choriocapillaris layer and the photoreceptor plateau was calculated for all methods in three normal subjects (Eq. (9)). The results are presented in Tables 5 and 6. In case of single volume imaging, a mean CNR was calculated from 4 single volume data sets acquired in each subject. The standard deviation of this mean was calculated as a measure of error. In case of the volume averaging one mean volume was generated from 4 single data sets in each subject. Thus, only a point estimate of the contrast to noise ratio was calculated. Fig. 12. Comparison of OCTA depth profiles obtained in single (not averaged) volumes in subject N1. First column: amplitude-decorrelation method. Second column: phase-variance OCT. Third column: speckle-variance OCT. First row: full spectrum, single volume data. Second row: split spectrum, single volume data. Horizontal, dashed lines indicate the level of the "photoreceptor plateau." The Gaussian curves indicate depth windowing of the data applied for calculation of the CNRs reported in Table 5. Blue curve: "photoreceptor plateau." Purple curve: choriocapillaris. The PV OCTA method provided higher CNR than the two other methods in the single volume imaging results, as well as in the volume averaging method. Application of the splitspectrum technique tends to increase the CNR for each of the methods and within each subject. However, the difference is smaller than the error. Increased CNR value larger than single volume errors is achieved within the techniques (full spectrum and split spectrum) and within the subjects by averaging multiple volumes. Figures 14 and 15 compare choroidal images obtained with the three OCTA methods. The images of choriocapillaris obtained from single data sets have a "granular" appearance, only vaguely suggesting a mesh-like structure. Volume averaging yields a structure more convincingly resembling the appearance of choriocapillaris as imaged in histology samples [3]. A comparison between the three methods also shows differences in the influence of the retinal blood flow on the choriocapillaris imaging. In the PV OCTA method, the retinal vessels provide strong OCTA signal ("Doppler traces") while the SSADA and SV OCT method do not show the retinal vessels in images obtained from Subject N1.   Table 6. Blue curve: "photoreceptor plateau." Purple curve: choriocapillaris. The images of Sattler's layer (Fig. 15) show similar trends, i.e., volume averaging improves the visualization of vessels. The PV OCTA method reveals the most detailed image of this vascular layer. Most interestingly, in the OCTA projection obtained from the average volume ( Fig. 15(h)), Sattler's layer vessels are visualized with the choriocapillaris structure visible in the background (which may suggest the vessels supplying choriocapillary lobules). Figure 16 shows results of split spectrum PV OCT imaging with application of the volume averaging in three normal subjects. The comparison of the OCTA depth profiles shows the presence of a sharp choriocapillary peak in all subjects. However, subject N1 has significantly stronger signal originating from the choroid than subjects N2 and N3. The locations of the Gaussian windows applied for generation of the OCTA projections show consistency in selection of the choriocapillaris location which coincide with the sharp peak. However, the best visualization of Sattler's layer was achieved at different depths for all three subjects. In the oldest subject (N1, age 63) the Gaussian window had to be placed closer to the choriocapillaris peak (15µm) than in the two younger subjects (N2 age 35: 29 µm and N3 age 29: 39µm).

Between-subject comparison of PV OCT imaging of the choroid
In the OCTA projections, the choriocapillaris is clearly visualized in all cases similar features. The mesh-like structure is apparent in all images. More variability is visible between the subjects in the images of Sattler's layer. Subject N2 has the poorest visualization of the vessels which have a very fragmented appearance. Images of subjects N1 and N3 show clearer vessel structure.

Discussion and conclusions
The OCTA methods have the capability to visualize all vascular layers present in the retina. The quality of visualization, however, is dependent on the scattering properties of the anatomical layers. The sharpest image of the vascular system was obtained at the edge between the inner nuclear layer and the outer plexiform layer. Due to high transparency of the tissue, the outer capillary plexus was clearly visualized even in the intensity projections. First column: amplitude-decorrelation method. Second column: phase-variance OCT. Third column: speckle-variance OCT. First row: full spectrum, single volume data. Second row: split spectrum, single volume data. Third row: average of 4 volumes obtained with the splitspectrum data. Image location: 6° nasal, 4° inferior from the fovea. Image sizes: 1.8 x 1.8 mm.
In highly scattering tissues, however, the OCT angiography methods show superior performance over the standard intensity imaging. OCTA techniques provide visualization of the inner capillary plexus and the larger retinal vessels with clarity not attainable by intensity imaging. The most challenging imaging conditions were found in the nerve fiber layer, in which only fragmented images of capillaries were generated.
In the visual comparison between the images obtained with the three implemented OCTA techniques, the amplitude-decorrelation method rendered the smoothest images of the retinal vasculature (less grainy appearance than the speckle-variance and phase-variance methods). Application of the split spectrum approach improves the performance of all the methods. However, a major improvement in the continuity of vascular images was obtained with the volume averaging.
Imaging of the choroid is challenging due to its location underneath highly scattering RPE. Yet, the possibility of visualizing its vascular layers may be clinically important as it plays a critical role in maintaining proper functioning of the RPE and photoreceptors. Overall, OCTA imaging of the choroid demonstrates that the OCTA methods can provide images of the choriocapillaris and Sattler's layer. However, vessels located in Haller's layer and below do not produce OCT signal. Therefore, at present, they can be visualized only in the shadowgrams (the lack of signal acts as a contrast in the intensity imaging) and are sometimes displayed in inverse intensity gray scale.
In OCTA imaging of the choriocapillaris, careful data processing is required. The thickness of this layer measures from 5 µm to 18 µm at the posterior pole of the eye [1,80]. To avoid sectioning through the choriocapillaris layer, the visualization in en face projections should closely follow the shape of the outer edge of the RPE (Bruch's membrane) in the healthy eye. In OCT angiography imaging the choriocapillaris layer produces several characteristic features. In the OCTA B-scans, a bright band of flow signal abuts Bruch's membrane. In average OCTA depth profiles, the choriocapillaris gives rise to a very sharp edge clearly marking the boundary between Bruch's membrane and the choroid. In OCTA en face projections, a characteristic meshwork of flow signal is visible resembling Imaging of Sattler's layer seems to vary more among individual eyes than imaging of the choriocapillaris. In two of the three imaged normal subjects vessels were clearly visualized; in the third, the appearance was fragmented. The position at which the clearest image of Sattler's layer vessels was generated in en face projections differs between the subjects. In subject N1 (age 63) the projection was located 15 µm below the choriocapillaris peak. In subject N2 (age 35) the position was 29 µm, and in subject N3 (age 29) the distance from choriocapillaris was 30 µm. Although this may suggest age-related differences, we do not have enough experimental evidence to evaluate this hypothesis. Also, imaging of this layer shows differences in the performance of the compared OCTA methods. In our implementation, the PV OCTA method rendered the clearest images of Sattler's layer vessels.
Wide-field intensity imaging of Haller's layer provided visualization of large choroidal vessels. The vessel contrast depends however on the light penetration into the choroid in individual eyes which also depends on the condition of ocular tissues. In diseases compromising the function of RPE and choriocapillaris, increased depth penetration of the imaging beam allows for better visualization of deep choroidal vessels [38,79,90], as demonstrated in the comparison between normal, AMD and GA cases. In projection images generated at the scleral level we also observed vessels located below Haller's layer. Previous reports exist in which imaging artifacts were identified in the OCT imaging of sclera [14]. However, the orientation of the features in our imaging and their connection with Haller's layer vessels are inconsistent with the artifacts and suggest vessels, possibly short ciliary arteries or veins of lymphatics [82,87,91]. In general, careful analysis of features present in images of tissue located underneath highly scattering layers is required to avoid erroneous interpretation of artifacts generated by the imaging methods as real structures.
The comparison between the amplitude-decorrelation, speckle-variance and phasevariance OCT angiography methods was performed by analysis of the average depth profiles and by calculation of the CNR in the imaging of choriocapillaris. In our implementation, the three methods show differences in the suppression of the static tissue noise. The strength of the OCTA signal originating from the photoreceptor layer measured relative to the height of the choriocapillaris peak served as an estimation of performance. The weakest OCTA signal from photoreceptor IS/OS junction was observed in the phase-variance method suggesting the most effective suppression of static tissue noise, even with "Doppler traces/shadows" giving more contribution to the signal intensity in the RPE complex than the other methods. Also, in calculations of the CNR for choriocapillaris, the phase-variance method provided the highest values.
Application of the split-spectrum method improves the CNR in all methods. However, in our implementation, volume averaging surpasses the improved performance of the splitspectrum technique in all three compared OCTA methods. Poorer results of application of the split spectrum method may be caused by non-optimized splitting criteria. In general, higher separation of the centroids of the spectral bands and lesser overlap between them should provide better results [92]. Our conservative choice (wide and overlapping bands), may not be sufficient to observe significant improvement in the OCTA images. In the papers reporting optimization of the split-spectrum procedure [48] a number of splits N = 11 were found to optimize certain metrics (signal-to-noise ratio and vessel connectivity) used to compare visualization of vessels. However, the bandwidth narrowing factor was 0.23 of the FWHM of the original spectrum. If we applied this result, we would decrease our axial imaging resolution to ~26 µm compromising the main advantage of the OCT technique -depth resolution.
Although only three healthy subjects were tested, the results were consistent among them, and provide insight as to what kind of differences may be expected from a larger sample with the compared OCTA techniques. For example, techniques that perform well in imaging of the retinal vasculature may not perform equally well in the imaging of the choroid, or may need additional optimizations of the imaging or data processing methods. Similarly, methods showing higher sensitivity in detection of weak flow signals from the choroid may introduce stronger flow induced artifacts in the images of the retinal vasculature obscuring the visualization of vessels ("Doppler traces/shadows" in scattering tissue underlying vessels).
A multitude of OCT angiography methods has been developed in recent years generating questions regarding their advantages and limitations in ophthalmic research and clinical diagnostics. In this paper we provided an overview and comparison of three commonly used techniques. Although differences exist among them, no single universal method can be identified, yet. Each method may have advantages in imaging of specific tissues and therefore their usage may depend on specific applications. It is clear however that volume averaging aided by motion correction methods will be critical for improved clinical OCTA systems for comprehensive chorio-retinal vascular mapping in health and disease.

Appendix 1. Implementation of the amplitude-decorrelation and speckle-variance OCTA methods
In amplitude-decorrelation OCTA and speckle-variance OCTA methods the amplitude (linear scale) of the complex OCT signal is used to visualize flow. In our implementation, the common processing steps used for both methods are as follows.

Thresholding of the intensity images
All OCTA methods use thresholding of the intensity data to ignore the noise pixels that otherwise obscure the information about flow in the final visualizations (noise gives much higher values of decorrelation and variances than the areas where flow occurs). Mean intensity and standard deviation was calculated for the entire volume data set. Voxels with intensity lower than 50% of the mean to standard deviation ratio were replaced with zero values. This threshold was found empirically by visual inspection of the resultant OCT angiography images and was applied to remove intensity noise which would introduce erroneously high variance and high decorrelation values in further processing steps.

Registration of repeated MB-scans
Small motions of the eye and mechanical instabilities of the OCT apparatus cause shifts of the intensity distribution in the acquired cross-sections within the MB-scan series. This introduces erroneously high OCT angiography signal to the images of static tissue (OCTA noise) and obscures the flow signal from the vessels. Significant motion which entirely changes the intensity distribution of the speckle field (e.g. saccadic motion of the eye) cannot be corrected by simple cross-correlation. Artifacts caused by the shift of the speckle field without significant changes in the intensity distribution can be reduced by correlation of the B-scans. The correlation was performed prior to implementation of the OCT angiography computations to reduce the loss of the flow signal.

Calculation of the OCTA metric
The data processing steps which differ for the speckle-variance and amplitude-decorrelation methods are dictated by the calculations required by the metrics inherent to these methods. OCT angiography images were generated using Eq. (4) for SV OCTA and Eq. (7) for AD OCTA methods. Care needs to be taken while computing the amplitude decorrelation values to avoid dividing any pixel by a zero value.

Removal of the noisiest components
B-scans introducing the highest noise were removed from the MB-scan series. In the specklevariance method, "trial" angiography images were calculated from groups of 4 B-scans out of 5 acquired within the MB-scan group, by removing 1st, then 2nd, 3rd, 4th and 5th B-scans. Out of 5 "trial" angiograms the one with the lowest mean variance signal was selected as the final output. This procedure was repeated for each MB-scan group in the acquired volume.
In the amplitude-decorrelation method, pairs of B-scans were used to calculate the OCT angiography images. For each group of 5 B-scans within the MB series, 4 pairs of B-scans (1-2, 2-3, 3-4, 4-5) were used to compute the decorrelation images. The pair that produced the highest mean decorrelation signal was removed as the noisiest component prior to averaging of decorrelation images according to Eq. (7). which the eye can move involuntarily introducing motion artifacts to the OCT angiography images. Motion which does not cause full decorrelation of the OCTA B-scans can be corrected by using the method described in Appendix 1, section 2. Large and rapid motions of the eye (e.g. saccadic motion) introduce uncorrectable motion artifacts visible in OCT en face projections as horizontal white lines ( Fig. 17(b)) which correspond to angiography B-scans with motion erroneously detected in the entire image of the retina (Fig. 17(e)). Bulk motion-affected OCT angiography B-scans were removed from the data set by automatic detection of cross-sections with the mean intensity exceeding an empirically selected threshold. Mean intensity was calculated for every cross-sectional OCT angiography image in the acquired data volume resulting in a vector of mean values. An absolute value of this vector derivative was computed showing the abrupt changes in the mean intensity as "spikes" (Fig. 17(a)). OCT angiography B-scans corresponding to "spikes" with magnitude larger than standard deviation of the mean were removed as motion affected outliers. If the distance between the "spikes" was less than 5 pixels, the entire block of B-scans between the peaks was removed as motion affected. This procedure broke the volume data sets into blocks of saccadic motion artifact-free B-scans. The location coordinates of these blocks within the data volume were stored in a separate vector. Slower, sheer, torsional, rotational, etc. motion may be present in the blocks causing distortion of the imaged vasculature but without significantly obscuring the OCT angiography signal. Correcting for these types of motion is outside the scope of this study.

Flattening of the retinal images
The arrangement of choroidal layers in the normal eye follows the curvature of the outer retinal complex (photoreceptors layer and retinal pigment epithelium). Flattening to this layer was performed to visualize the choroidal vasculature. The flow chart of the implemented procedures is shown in Fig. 18.  a. Median filtering was applied in each B-scan to remove noise and smooth out the intensity distribution (reduce speckle pattern). The median window size was 25 pixels laterally and 5 pixels axially providing high smoothing along the retinal layers.
b. Intensity thresholding was implemented to find the boundary between the outer nuclear layer (highly transparent in OCT imaging) and the inner/outer photoreceptor segments junction (IS/OS junction; highly scattering) which reflects the shape of the outer retinal complex. In each A-scan the two strongest intensity thresholds usually corresponded to the nerve fiber layer (lower position index) and the IS/OS junction (higher position index). The second one was selected as the IS/OS junction location estimate.
c. To validate the estimated locations of the IS/OS junction the mean of all found positions was calculated for each B-scan. The data point closest to the mean served as a starting reference point for validation of the remaining points. If the next position estimate was within a range (−15; 15) pixels relative to the reference point, it was considered as valid and served as a new reference against which the next position estimate was checked. Otherwise, the estimated position was rejected as not valid and the last valid point remained a reference for the next estimate. The same procedure was applied to position estimates preceding the starting mean data point.
d. The validation procedure was performed in each B-scan in the acquired volume. The output of the IS/OS junction position search was a matrix of the edge location coordinates (X, Y, Z) which was used to perform a surface fit. To match the curvature of the outer retinal complex, Zernike polynomials up to the 9th order were used in the 3D surface fitting procedure. The obtained surfaces were used to flatten the intensity and OCT angiography images.

Transverse motion correction and volume averaging
Four 3D data sets were used in transverse motion correction and volume averaging. The flow chart of implemented procedures is shown in Fig. 19. b. Cross-correlation of overlapping en face projections of motion-free data segments was performed to correct for transverse eye motion. The largest segment in all used data sets was selected to initialize the motion correction procedure. The next largest en face segment with at least 20% overlap was correlated with the initial one. The 20% overlap condition was checked based on the original location of the segment in its corresponding volume. The overlap condition implied that the two segments originated from different volumetric data sets. If the automatic correlation procedure was erroneous, the user had an option to manually correct the position of the segment. The corrected transverse coordinates (X, Y) of the segment were stored in the remapping matrix. The two overlapped en face segments were then averaged (the summed intensity was divided by 2 in the area of intersection, and by 1 in the area of symmetric difference) and served as a template for correlation with the next largest segment. As the procedure progressed, subsequent segments were merged with the expanding template until the pool of all segments was exhausted. The output was the remapping matrix which was used to shift all the 3D data blocks to the motion corrected XY locations.
c. The XY remapped data blocks were registered to remove axial displacements between the original volumetric data sets. B-scans from the centers of overlapping blocks were correlated and the blocks were shifted accordingly.
d. In the final step, the XYZ registered data blocks were averaged providing motion corrected volumes.