Gabor optical coherence tomographic angiography (GOCTA) (Part I): human retinal imaging in vivo

: Recently, parallel high A-line speed and wide field imaging for optical coherence tomography angiography (OCTA) has become more prevalent, resulting in a dramatic increase of data quantity which poses a challenge for real time imaging even for GPU in data processing. In this manuscript, we propose a new OCTA processing technique, Gabor optical coherence tomographic angiography (GOCTA), for label-free human retinal angiography imaging. In spectral domain optical coherence tomography (SDOCT), k-space resampling and Fourier transform (FFT) are required for the entire data set of interference fringes to calculate blood flow information in previous OCTA algorithms, which are computationally intensive. As adults' eye anterior-posterior radii are nearly constant, only 3 A-scan lines need to be processed to obtain the gross orientation of the retina by using a sphere model. Subsequently, the en face microvascular images can be obtained by using the GOCTA algorithm from interference fringes directly without the steps of k-space resampling, numerical dispersion compensation, FFT, and maximum (mean) projection, resulting in a significant improvement of the data processing speed by 4 to 20 times faster than the existing methods. GOCTA is potentially suitable for SDOCT systems in en face preview applications requiring real-time microvascular imaging.


Introduction
Optical coherence tomography (OCT) [1] technique, proposed in 1990s, is an emerging imaging modality for medical diagnostics and treatments. Due to the advantages of non-invasiveness, high resolution and high imaging speed, OCT is widely used for various tissue, e.g. human retina, brain, cardiology and dermatology. In addition to microstructural imaging, OCT based microvascular imaging algorithms are also used widely in medical imaging and play an increasingly important role. The first algorithm for extracting blood flow information is optical Doppler tomography (ODT) [2,3] or color Doppler OCT (CDOCT) [4] which is able to calculate the axial velocity component of moving scattering particles.
Morphological OCT microvasculature imaging, collectively termed OCT angiography (OCTA), was also developed. In general, OCTA algorithms available now can be divided into two categories according to processing mode. The first is inter-line mode, such as Doppler variance phase resolved (DVPR) [3,5], intensity-based modified Doppler variance (IBDV) [6,7], optical micro-angiography (OMAG) [8]. For inter-line mode, the blood flow information was extracted from one frame of interference fringes at each position. For DVPR and IBDV, the statistical information of a small window was calculated to contrast microvasculature, which needs high A-line density. For OMAG, a piezo-stage was used in reference arm for modulating of interference fringes, which increased complexity of OCT setup.
The second processing mode is inter-frame, which extracts blood flow information from multi-frames of structural images at each position, such as phase variance OCT (PVOCT) [9][10][11], speckle variance OCT (SVOCT) [12][13][14][15], correlation mapping OCT (cmOCT) [16][17][18], split-spectrum amplitude-decorrelation angiography (SSADA) [19] and differential standard deviation of log-scale intensity (DSDLI) [20], and ultrahigh sensitivity optical microangiography (UHS-OMAG) [21][22][23][24]. For this mode, the sensitivity for microvasculature detection can be improved due to the time interval between two frames is longer than that between two A-scans, but the motion artifacts is also more significant due to increase of time interval. PVOCT, SVOCT, cmOCT, SSADA, and DSDLI obtain blood vessel contrast by calculating statistical information from either phase or intensity images in spatial domain. PVOCT calculates the variance of phase difference between two frames. SVOCT and DSDLI calculate the variances of intensity and the differential intensity between two frames, respectively. Both cmOCT and SSADA calculate the decorrelation coefficients, but in SSADA, the full spectrum is divided into four sub-bands to improve microvascular image quality. For UHS-OMAG, OMAG algorithm is performed in the slow scanning direction and blood flow signal is calculated from both amplitude and phase signals, resulting in an improvement of sensitivity.
Recently, parallel imaging [25][26][27] and wide field imaging [28] have become more prevalent, resulting in a dramatic increase of data quantity which poses a challenge for real time imaging even when using GPU for data processing. For all of the above mentioned algorithms, the blood flow information is obtained from spatial domain. SDOCT systems all require k-space resampling, dispersion compensation, Fourier transform (FFT), and maximum (or mean) projection (MIP) for reconstructing en face microvascular images, which appears to be the most useful display mode for clinical applications. Some of these steps require long processing times, which poses challenge for real-time imaging even when using GPUs for data processing. For clinical applications such as retinal imaging, we recognize that OCTA images are typically used as en face image sets for clinical decision making, such as identifying an area of microvascular abnormality, after which depth resolved information, such as cross-sectional structural OCT images of the retina at the particular region, are reviewed. Therefore, rapid en face OCTA image display, at the time of scanning, may be advantageous to screen retinal pathology as well as to focus detailed examination to a smaller region of interest. Such capability may also be useful for less cooperative patients where motion artefacts degrade OCTA images. In such scenarios, rapid en face OCTA may allow immediate feedback and re-scanning.
It is thus interesting to reflect that most existing OCTA algorithms go through the computationally intensive steps of depth resolved image processing, and in the last steps, perform intensity projection in the depth direction and remove the depth information. In this work, we proposed a Gabor optical coherence tomographic angiography (GOCTA) algorithm extracting blood flow information from interference fringes directly without the time consuming steps mentioned above, which can decrease data processing time substantially.  Figure 1 illustrated the data processing steps of GOCTA. In SDOCT, the interference fringes between the light backscattered from sample and reflected by reference mirror are detected by spectrometer camera. The three dimensional (3D) data set of spectral fringes can then be acquired by scanning of the x-and y-galvo-mirrors. The direct component (DC) of the interference (auto-correlation of reference beam) can be measured by blocking sample arm and the auto-correlation of sample beam is negligible. After subtracting the DC component, the captured signal can be simplified by

Method
where x and y represent the scanning directions of the two galvos, λ is wavelength, ( ) S λ is the power spectral density of light source, ( , , ) x z y R and r R are the backscattering coefficient of sample and the reflectivity reference mirror, respectively. s γ and r γ are the input power in the sample and reference arms, n is the refractive index, z represents depth, ( ) are the initial phase and the dispersion mismatch between sample arm and reference arm. In the case of moving particles, the amplitude and the frequency of the fringes vary with time.
For two consecutive B-scans acquired from the same position, the amplitude or frequency of the components corresponding to moving particles is different. Subtracting the two B-scans, the components corresponding to static tissue can be removed, and the resultant signal originates from the moving particles. The differential signal can be expressed by As shown in Fig. 2, human eye as an optical system has a curved image plane on the retina near the fovea, with the optical nerve head in the vicinity. The anterio-posterior (AP) diameter of an emmetropic human adult eye is approximately 22 to 24 mm, which is relatively invariant between sex and age groups [29]. While the transverse diameter varies more according to the width of the orbit, for the area near the fovea and optical nerve head, the curvature can be approximated by the AP diameter. For the work presented below, we have chosen an AP diameter of 22 mm (or radius of curvature of 11 mm) for the human eyes. In GOCTA, the orientation of retina is needed for generating Gabor filters and evaluated by the spherical model which can be expressed by where (x 0 , y 0 , z 0 ) and R are the center position and the radius, respectively, z s is the depth of retinal surface in structural images. To calculate the center location (x 0 , y 0 , z 0 ), three sets of (x, y, z s ) are needed. Here, we performed FFT on the three A-scans at corners and calculated the depth of surface. The retinal surface z s (x, y) can be then obtained by Eq. (3), approximating the region marked by the dashed box in Fig. 1. Within the measured interference fringes, the sample information at different depth is modulated by different frequency component. As Gabor filter is a linear filter, the frequency component within a specific frequency range can be obtained directly by convolution [30], which is equivalent to multiplying a Gaussian function in spatial domain. For example, Δ are the depth and depth range respectively. Taking the refractive index and round optical path into account, the filter can be obtained by performing FFT on the above mentioned Gaussian function and expressed by where M is the pixel number of CCD, '' ( , ) mean I x y is the mean value of each A-scan of the filtered fringe. By calculating GOCTA signal for each position in the 3D data set of spectral fringes, en face microvascular images can be directly obtained.

OCT imaging system
We performed GOCTA algorithm on the data set of healthy human eye from a commercial SDOCT (AngioVue, OptoVue Inc.) to verify its performance. This system operated at a center wavelength of 840 nm with the axial resolution and lateral resolution of ~5 μm and ~15 μm, respectively. The A-scan rate is 70,000 A-scans per second. In this work, the scanning range is 3 × 3 mm and each position was scanned twice.

Human retina imaging
We performed retinal OCT scanning on ten healthy volunteers. Example data for two local regions (optical nerve head region and fovea region) are shown in Figs. 3 and 4, respectively. The scanning ranges were 3 × 3 mm with 608 × 304 A-scans. SVOCT, UHS-OMAG and SVOCT algorithms were performed on the same data set to calculate microvascular images for comparison and the en face images were obtained by using mean projection within the depth range, same as the result obtained by Gabor filters. All of the en face microvascular images were calculated within depth of 0 -350 μm below the retinal surface. Signal to noise ratio (SNR) and contrast to noise ratio (CNR) of the en face micro-vascular images were also calculated for quantitative comparison, SNR and CNR were calculated by To quantitatively assess the microvascular and background signals for comparison, we double threshold the marked regions to obtain the masks for dynamic signals (red) and background (blue), as shown in Fig. 3(c) and Fig. 4(c), which include vessels of different sizes. The results demonstrate that GOCTA can provide comparable image quality compared to the other three algorithms in the vicinity of both the optical nerve head and fovea regions, as shown by the comparable SNRs and CNRs. The output images from commercial SSADA algorithm use duplicated optical scanning in two directions (x and y) with post-processing applied on these two data sets to suppress motion artifact. We measured and calculated SNRs and CNRs on all ten volunteers' data sets from healthy volunteers and the average and standard deviation of SNRs and CNRs for GOCTA, SVOCT, UHS-OMAG and SSADA are 25 ± 2, 23 ± 2, 23 ± 2, 22 ± 1 and 14 ± 2, 9 ± 2, 11 ± 3, 9 ± 2, respectively. The results demonstrate that the proposed GOCTA algorithm can provide comparable image quality.

Comparison of data processing time for the four algorithms
It is important to recognize that the commercial system output of SSADA images as shown in Fig. 3(b) and 4(b) are based on Fig. 3(f) and 4(f), upon which corrections of motion artefact, sub-regional registration, and additional contrast enhancement are applied to two directions (x and y) of scanning data set. We noted the GOCTA images quality, both qualitatively and quantitatively, were at least on par (or possibly superior) in comparison with the unidirectional scan results of SSADA. These subsequent data processing techniques on the commercial system are proprietary and thus not compared for computational complexity and time. In principle, these would be equivalent when applied to any en face data set.
The main advantage of our newly proposed GOCTA algorithm is the speed of processing. We processed the data set on the same computer using published SVOCT, UHS-OMAG, and SSADA algorithms, in MatLab. Note that in order to obtain data sets used to post-process the commercial SSADA image, both scanning in the x and y directions were performed and the SSADA algorithm must be repeated, which doubled the numerical processing time. The data processing was accomplished on a laptop (CPU: i7-4720HQ, memory: 16 G, GPU: NVIDIA Geforce (GTX 970M), operating system: windows 8.1). The data processing time for each 2 B-scans from the same position was calculated and the results are shown in Table 1. GOCTA is the only algorithm that can directly provide en face microvascular images without FFT. In SVOCT, UHS-OMAG, and SSADA, the steps of k-space resampling, dispersion compensation and FFT are computationally costly, resulting in 6, 4, and 20 times slower than GOCTA, respectively. Since GOCTA does not require resampling, dispersion compensation, and FFT, the total processing time decreases dramatically. Using GPU based parallel computing library on Matlab, we also measured the data processing time for the entire 3D (608 × 2048 × 304) data set and the results are shown in Table 2. Note these results are only illustrative as additional GPU algorithm optimization is possible outside of the Matlab programming environment, however, the comparison here serves to show the different components of the overall computational can be GPU accelerated to different extents. We simulated the data processing on GPU, assuming the acquired interference fringes were transferred to GPU frame by frame (2048 × 304), as each B-mode image data set became available. For such scheme to work, the real-time OCT imaging system will require the three corner A-scans to be performed first, obtain the necessary spherical approximation of the retina, and apply the GOCTA as well the other processing algorithms for comparison. The processing time for these three A-scans is negligible with respect to the overall processing time.
It is also important to note that the steps of k-space resampling, numerical dispersion compensation, and image alignment required both matrix amplitude and matrix indices adjustment using algorithms such as spline fitting, no GPU based parallel computing Matlab library was readily available. Hence these were kept as CPU operations in the current analysis but could be further improved outside of the Matlab environment. Nevertheless, since the overall computational complexity of GOCTA was simpler than SVOCT, OMAG, and SSADA, the above analysis showed GOCTA was faster to compute under GPU acceleration. In this work, k-space resampling was accomplished by using cubic spline interpolation.

Limitations and conclusion
In the proposed algorithm, the Gabor filter parameters were deliberately chosen such that large number of zeros was encountered, thus simplifying the computational complexity and reduce the time needed for the convolution in digital filtering. In this work, the microvascular images within depth range of 350 µm (10% of the total OCT ranging depth) below spherically fitted retinal surface are calculated for analyzing and comparison. In this case, the non-zero segment length of Gabor filter Eq. (4) is just 16 pixels, resulting in a substantial decrease of computation complexity.
Due to the proposed GOCTA algorithm skipping the step of calculating depth resolved structural images (in the z-direction), one of the limitations is that the structural image alignment in the z-direction cannot be performed for motion artefact removal. However, xand y-direction based en face image registration and alignment can still be applied (e.g., [31,32]). The curvature of lens system can affect the accuracy of the evaluated retinal orientation, for slight curvature, GOCTA's images will not be affected due to the depth range of Gabor filter is 10% of the total OCT ranging depth. In the case of significant curvature, the relative shifting distance at each pixel can be obtained by scanning a mirror and the evaluated retinal orientation can be compensated in software. Another limitation is that GOCTA can only provide en face images, but for preview images, as our results shown above, it can still provide useful diagnosis information for the ophthalmologist. Finally, we noted commercial systems may use algorithms other than cubic spline interpolation for k-space resampling for processing efficiency with trade-off in image quality.
As the results shown, the SNRs and CNRs obtained by GOCTA are slightly higher than the other three algorithms. The reason might be that the proposed algorithm uses a large range of frequency component (the sample information within depth range of 2 z z δ − Δ to 2 z z δ + Δ in spatial domain) to calculate the blood flow information, which is more robust compared to the other three algorithms where only the sample information at the same depth is used and then perform maximum (mean) projection to generate en face microvascular images.
In conclusion, we proposed a novel Gabor optical coherence tomographic angiography (GOCTA) algorithm to obtain the microvasculature on human retina on a standard ophthalmic SDOCT system. Obviating the need of resampling, dispersion compensation and FFT, the proposed algorithm can achieve a 6, 4 and 20 times of the data processing speed compared to SVOCT, UHS-OMAG and SSADA, respectively. GOCTA is ideally suited for SDOCT systems used in wide field scanning, ultra-high spectral resolution or parallel high A-line speed applications where large data amount is generated. For real time imaging, GOCTA can also be performed on graphics processing units (GPUs) to increase data processing speed further. We consider the proposed algorithm potentially useful for calculating the preview OCTA images as the first line en face display for the clinician, and improve the efficiency of disease screening and diagnosis in busy clinical environment.