Regression-based algorithm for bulk motion subtraction in optical coherence tomography angiography

: We developed an algorithm to remove decorrelation noise due to bulk motion in optical coherence tomography angiography (OCTA) of the posterior eye. In this algorithm, OCTA B-frames were divided into segments within which the bulk motion velocity could be assumed to be constant. This velocity was recovered using linear regression of decorrelation versus the logarithm of reflectance in axial lines (A-lines) identified as bulk tissue by percentile analysis. The fitting parameters were used to calculate a reflectance-adjusted upper bound threshold for bulk motion decorrelation. Below this threshold, voxels are identified as non-flow tissue, their flow values are set to zeros. Above this threshold, the voxels are identified as flow voxels and bulk motion velocity is subtracted from each using a nonlinear decorrelation-velocity relationship previously established in laboratory flow phantoms. Compared to the simpler median-subtraction method, the regression-based bulk motion subtraction improved angiogram signal-to-noise ratio, contrast, vessel density repeatability, and bulk motion noise cleanup in the foveal avascular zone, while preserving the connectivity of the vascular networks in the angiogram.

Software and hardware solutions exist to remove the effects of eye motion on OCTA. Microsaccadic motion can be corrected by registration methods [23][24][25][26][27] that retrieve the signal missing in corrupted frames from other scans. In addition, tracking-assisted scanning schemes [28][29][30][31][32] have been implemented in commercial systems to prevent data recording during a microsaccade. Also, inter-volume optical microangiography [33] operates at an ultrahigh volumetric imaging rate in order to detect flow by recognizing signal variations between adjacent volumes instead of consecutive frames, showing potential to reduce the prevalence of motion artifacts. Yet, small-amplitude bulk motion caused by ocular drift, pulsation, tremors, or OCT system mechanical instabilities is currently undetectable by tracking systems and only partially attenuated by averaging registered images.
One way to reduce the effect of bulk motion is to co-register the reflectance images at each scanning position before computing the flow image. Methods that iteratively maximize cross-correlation of successive B-scans prior to OCTA processing have been used to estimate displacement and compensate bulk image shifts as well as global phase variations in the axial and lateral directions [34][35][36]. However, the actual three-dimensional nature of eye motion during scanning challenges this approach, which is limited to in-plane shifts. Another approach to solve the bulk motion problem is to approximate its velocity to the spatial average (along A-scan and B-scan directions) of the reflectance OCT signal variation between consecutive scans [37]. This solution requires acquisition of many frames of the same B-scan and assumes a low capillary density on the imaged sample.
In our split-spectrum amplitude decorrelation angiography (SSADA) algorithm, a simple bulk motion subtraction method was introduced to minimize the noise caused by motion between two consecutive B-scans. The decorrelation caused by bulk motion is estimated by the median decorrelation value within retinal tissue region, and then subtracted from all pixels at each B-frame of SSADA [9]. Although the median subtraction method improves the image quality, it constitutes an inaccurate approximation of the actual bulk motion contribution to flow signal.
In this paper, we propose a more accurate mechanism to remove the noise introduced by bulk motion. Using regression analysis; we first define the relationship between background bulk motion decorrelation and the logarithm of the OCT reflectance in a group of A-lines devoid of flow voxels within a segment of a B-frame. Then, we assign zero decorrelation values to the non-flow pixels in each segment where the respective voxels have values below a certain threshold derived from the relation between bulk motion and reflectance. A nonlinear model that relates decorrelation and flow velocity [38] is used to subtract estimated bulk-motion velocity from the measured velocity inflow voxels. Compared to the median subtraction algorithm, we demonstrate reduced bulk motion at avascular regions, improved image quality and vessel density repeatability while preserving vascular connectivity.

Data acquisition
Volumetric OCTA scans of the macula and the optic disc of healthy subjects were acquired using a wide-field 200-kHz OCT system. The system relies on a swept-source configuration that utilizes a tunable laser (Axsun, Inc., Billerica, MA, USA) operating at 1044 nm center wavelength with a 104 nm tuning range. It has a lateral resolution of 12 µm and an axial resolution of 7.5 µm, which is increased to 22.5 µm on the flow image. Participants were recruited at the Casey Eye Institute of OHSU and informed consent was obtained. The protocol was approved by the Institutional Review Board/Ethics Committee of OHSU and adhered to the tenants of the Declaration of Helsinki in the treatment of human participants. The scan pattern consisted of 800 B-frames composed of 850 A-lines and located at 400 raster positions (y-priority), covering an 8x6mm 2 area. Each B-scan was completed in 4.75 ms.

Data processing
Data was processed using the Matlab 2013a release (Mathworks, Natick, MA, USA). Crosssectional OCT images were generated and B-frames acquired at the same raster position were averaged. OCTA data was calculated using the SSADA algorithm [1]. Segmentation of the vitreous and inner limiting membrane (ILM) interface, outer plexiform layer (OPL) and outer nuclear layer (ONL) interface as well as Bruch's membrane and choroid interface was automatically performed by a directional graph-search algorithm developed by Zhang et al. [39]. The region contained between ILM and OPL is the vascularized inner retina. En face images were generated by maximum projection of the flow signal within the inner retinal slab. An expert grader manually corrected segmentation errors.
Frames affected by microsaccadic artifacts were recognized on the en face OCTA image by simply identifying summed flow signal above a threshold [28], and excluded from analysis. Then, we recognized the relationship between bulk motion contribution to decorrelation signal and the local reflectance. Because decorrelation is not linearly related to velocity, the subtraction of the bulk motion contribution must be performed on the velocity domain rather than in the decorrelation domain. Non-microsaccadic frames were divided into 5 segments where local reflectance does not vary abruptly. The segments are small enough that the axial-lines within them could be acquired within a small enough time frame (less than 1 millisecond) that bulk motion could be considered to be approximately constant within the segments. Then, we developed a method to estimate the bulk motion velocity from frame segments.
At each segment, we selected a group of A-lines devoid of inner retinal flow voxels. For this purpose, A-lines crossing large vessels were first identified by a binary large vessel mask of the en face angiogram and excluded from the following analysis. The mask was constructed by successively applying amplitude thresholding to remove the majority of capillary pixels, a morphological opening (erosion followed by dilation) to cleanup dispersed pixels, a Gaussian convolution filter to prevent holes in the middle of large vessels and a final binarization step. From the remaining A-lines, the first 10 percentile with the lowest decorrelation signal were selected, assuming these are composed by non-flow voxels only. The threshold chosen is approximately the upper boundary of inner retinal vessel densities among the healthy subjects that participated in the study. The decorrelation and natural logarithm of the reflectance data selected in the segment were converted into vectors and sorted by increasing order of reflectance values. Then, the vectors were divided into 20 equally sized bins and the reflectance and decorrelation values in each bin were averaged, forming the new vectors  1D). Then, the voxels in the k th segment with decorrelation below the threshold defined in Eq. (1) were set to zero. In frames affected by microsaccadic artifacts no A-lines were selected, the regression analysis step was skipped and the decorrelation values of all voxels were also set to zero. Fitting parameters in segments with less than five A-lines selected for regression analysis were substituted by the parameters of the neighboring segment within the same B-frame.
(1) Vascular voxels remain after thresholding, and the contribution of bulk motion , BM blood D to their flow signal was obtained using the median reflectance of the flow voxels within the frame segment ( blood R in Fig. 2A). An estimate of bulk motion velocity was calculated by a nonlinear model that related decorrelation and velocity in laboratory blood flow phantoms [38] (Eq. (2).
The saturation value sat D was calculated as 0.95 times the 99th percentile of the flow signal in the en face angiogram. brownian D was estimated assuming that it maintains the same proportion to sat D as the one reported previously [38]. The value of the saturation velocity sat v was approximated utilizing the experimental data obtained in Ref [38], adjusting for the difference in inter-scan period between the Avanti RTVue-XR (Optovue, Inc. Fremont, CA, USA) and the wide-field system used here. Equation (2) ) is modeled as a function of logarithmic reflectance R based linear regression coefficients m k and n k from the regression analysis of the k th frame segment. Voxels with decorrelation value below threshold D TH set at the 97.5 percentile of bulk motion (assuming normal distribution) are classified as non-flow voxels (D set to zero). The bulk motion-induced decorrelation in blood D BM,blood is calculated using the reflectance of blood blood R obtained from the median reflectance of flow voxels (voxel with D above D TH ).
Bulk motion velocity v BM is obtained from D BM,blood using Eq. (2). (B) The bulk motion velocity is subtracted from a voxel with uncorrected decorrelation value D 0 and corresponding velocity v 0 that is high relative to the saturation velocity according to the nonlinear curve relating flow signal D to velocity. Therefore subtracting v BM from v 0 to obtained corrected velocity v 1 and corrected flow signal D 1 has little effect. (C) The bulk motion velocity is subtracted from a voxel with uncorrected decorrelation value D 0 and corresponding velocity v 0 that is lower than the saturation velocity. Here subtracting v BM from v 0 to obtained corrected velocity v 1 and corrected flow signal D 1 has a larger effect than the previous example.
A flow diagram of the algorithm is shown in Fig. 3.

Comparison to median subtraction algorithm
The regression-based bulk motion subtraction method was compared with an earlier method in which the median decorrelation value of the segmented retinal region was subtracted from all retinal voxels within the B-frame segments defined above. The threshold used to distinguish vascular from non-flow pixels in original and median-subtracted angiograms was obtained from the decorrelation at the foveal avascular zone (FAZ) of en face projections by 2.33 , where D and σ represent mean value and standard deviation.

Comparison to inter-B-scan registration algorithm
An alternative approach to bulk motion removal a posteriori is to compensate the bulk motion between consecutive B-scans prior to computing OCTA flow signal. The benefits and tradeoffs of this pre-processing step in improving the artifact removal efficacy by the regression-based bulk-motion subtraction algorithm were evaluated.
For inter-B-scan registration we applied in SSADA a method similar to the one proposed by Lee et al [34]. Briefly, OCT images were first up-sampled by a factor of two in lateral (x) and axial (z) dimensions. An iterative routine maximizes the cross-correlation of the OCT signal of the two B-scans acquired at each position by rigid displacements of the second Bscan in two directions. After recognizing the optimal shift, it was applied to the second OCT image of each the 11 pairs generated by spectrum splitting [40]. Finally, the OCT images were down-sampled and the amplitude decorrelation signal was computed.

Image quality assessment
Bulk motion subtraction efficiency was assessed by the percentage of average FAZ signal remaining within the segmented retinal slab after processing the original data.
Parafoveal signal-to-noise ratio (SNR) was calculated by Eq. (4). The parafoveal annulus is concentric with the fovea, has an outer diameter of 2.5 mm and an inner diameter of 0.6 mm.
Vessel density was defined for en face projection as the percentage of area occupied by vascular pixels within the parafovea. Its coefficient of variation was used to assess inter-scan vessel density repeatability.
The RMS image contrast was defined as the standard deviation of the en face flow image, expressed in Eq.
Preservation of the vascular integrity after bulk motion subtraction was assessed by the vascular connectivity, defined in a skeletonized version of the en face angiogram by the percentage of flow pixels contained in groups larger than five. Skeletonization of en face OCTA is the process of converting vessels into lines with a one-pixel width. It was performed by the function bwmorph included in Matlab's Image Processing Toolbox. This function relies on an algorithm that iteratively removes pixels on boundaries of objects recognized on a binary image until the objects remain unchanged [41]. Since en face retinal angiograms contain vasculature of different dimensions (Fig. 4A), skeletonization was performed separately for large vessels and capillaries. If vessels with different calibers were not skeletonized separately, inaccuracies in the recognition of large vessels sometimes might be manifest (Figs. 4B-C). The p-value of a paired sample t-test evaluated the statistical significance of the results.

Results
Optic disc and macula from eight healthy volunteers were imaged. Scans containing microsaccadic artifacts crossing the parafoveal annulus were processed by the algorithm but excluded from the quantitative analysis.
A-lines containing microsaccadic artifacts and large vessels were successfully recognized (Fig. 5). The background A-lines selected at each segment avoided those containing vascular information (Fig. 5C). Segments where the fitting parameters had to be obtained from adjacent segments due to the small number of background A-lines available never represented more than 3% of the total. Segments with no A-lines selected were located at the microsaccadic artifacts. Changes in bulk motion background signal could be observed by simple inspection of wide-field frames affected by moderate bulk motion on en face images (Fig. 6A). The fitting parameters varied significantly along a B-scan (Fig. 6B), showing that segment partitioning helped reduce the effect of variations in local reflectance, allowing a more accurate bulk motion removal. Background signal at FAZ and bulk motion artefactual frames were effectively removed by the thresholding step (Figs. 7A-C), resulting in improved contrast between capillaries and background. Then, the brighter appearance of vascular pixels contained at artefactual lines was corrected by the step of bulk motion velocity subtraction (Fig. 7D). Larger subtraction from vascular voxels occurred in (non-microsaccadic) artefactual frames and smaller subtraction occurred at the flow voxels with decorrelation values near saturation (Fig. 7E). Compared to the prior median subtraction algorithm (Table 1), the regression-based bulk motion subtraction algorithm removed a larger percentage of decorrelation noise from the FAZ, achieved a greater improvement in vessel density measurement repeatability, a better signal to noise ratio for flow detection and a better RMS contrast. Two methods preserved similar vascular continuity. No improvement of RMS contrast was observed between mediansubtracted and original angiograms (p > 0.05), but significant improvement was found after regression-based subtraction. The vessel density did not see significant reduction from the original to the regression-based motion subtracted angiograms (p = 0.0687). The medianbased algorithm running time was 4.6 seconds while the regression-based algorithm took 20.3 seconds on CPU. A qualitative comparison of the two bulk motion subtraction methods is shown for en face angiograms of the macula and optic disc in Fig. 8 and for one representative cross-sectional B-frame in Fig. 9.   8. Qualitative comparison between median subtraction and regression-based bulk motion subtraction in an optic disc scan (A1-A7) and a macular scan (B1-B7).A1, B1 Unprocessed en face images. A2, B2 Scans after subtracting the median of the frame's retinal region. A3, B3 are obtained by subtracting A1-A2 and B1-B2 respectively. A4, B4 Scans after regressionbased bulk motion subtraction. A5, B5 are obtained by subtracting A1-A4 and B1-B4 respectively. A6 and A7 are close-ups of the 3.8mm × 1.7mm sections enclosed by dashed lines on A2 and A4 respectively.B6 and B7 are close-ups of the 3.4mm × 2.5mm sections enclosed by dashed lines on B2 and B4 respectively. The pre-compensation performed by inter-B-scan registration could help improve the artifact removal efficiency of the bulk motion subtraction algorithm ( Fig. 10(B) vs. Figure  10(D)). Addition of this step represented an average increase of 9 second of raw data processing per frame on CPU. Bright artifacts were reduced but not completely removed by inter-B-scan registration alone (Fig. 10C) and a background threshold was still necessary to identify vascular pixels. The pre-compensation could not remove the vessel density dependency on local reflectance by simple thresholding, while the regression-based thresholding proposed in this manuscript provided a correct spatial distribution of vascular pixels ( Fig. 10(G) vs. Fig. 10(H)). Additionally, the background-capillary contrast of the resulting image was slightly degraded compared to the case where no displacements are implemented prior to SSADA (Figs. 10(E) vs. Fig. 10(F)).

Discussion
Currently, correction of bulk motion contribution to flow signal in SSADA is performed by a median subtraction algorithm that underestimates the real bulk motion background and ignores the saturation of decorrelation at vascular voxels. In this manuscript, we propose a new method to increase the accuracy of bulk motion subtraction, improving image quality and enhancing the repeatability of vessel density quantification while preserving vascular integrity. This could improve image interpretation by reducing bright line artifacts, and make quantification of vessel density more accurate. Moreover, unlike in the median subtraction algorithm, vascular recognition from background does not require referencing to an avascular region such as the FAZ, allowing computation of the vessel density on scans outside the macula such as the optic nerve head and anterior segment. Although this algorithm is demonstrated for SSADA, it can be potentially applied to other OCTA implementations such as speckle variance [20] and optical microangiography [3].
The algorithm is based on three main steps: (1) regression analysis of the background decorrelation vs. reflectance curves, (2) reflectance-dependent thresholding and (3) subtracting bulk motion velocity from the vascular voxels to retrieve bulk motion-free decorrelation values. Both the thresholding step and the estimation of BM v rely on the accuracy of the linear regression routine, which assumes only background A-lines are selected. If voxels containing real flow signal were mistakenly chosen, they would contribute to larger flow signal at large reflectance voxels, overestimating the slope of the fitting curve. This would cause larger subtracted values, underestimation of the vessel density and degradation of the vascular integrity. We chose the first 10 percentile of A-lines with lower flow signal after removing the large vessel mask, given that the mean vessel density of macular scans did not exceed 90% for healthy subjects. This assumption underestimates the number of flow voxels. In optic disc scans, where some regions have a larger vessel density, the first five percentile of A-lines contained in the whole B-frame were selected. Averaging of voxel values in bins before regression analysis helped to minimize inaccuracies caused by occasional outliers. Implementation of inter-B-scan registration prior to SSADA alone could not correct the disadvantages inherent to fixed thresholding. On the other hand, implementation of the regression-based algorithm alone was more successful in recognizing vasculature from background and did so in a shorter time. If inter B-scan registration could be optimized to preserve or improve capillary contrast, it can be a supplementary step to further improve the BM removal efficacy.
This algorithm uses information available in a single scan to remove bulk motion background and small-amplitude artifacts. However, given the high prevalence of microsaccadic artifacts in OCTA, it would not eliminate the need for registration of more scans [23] in a clinical scenario. Although commercial systems have adopted various forms of real-time tracking to prevent the recording of microsaccadic motion [28,32], drifts between frames before and after microsaccades remain uncorrected, still necessitating registration of at least two volumes for truly artifact-free OCTA.

Conclusion
In summary, we have demonstrated a method that accurately subtracts bulk motion contribution to decorrelation of background and vascular voxels without affecting vascular integrity. The method recognizes the bulk motion decorrelation dependence on reflectance signal for background voxels; filters out background by a reflectance-dependent thresholding step, subtracts a bulk motion velocity value calculated by a nonlinear model that relates decorrelation and velocity, and finally retrieves the bulk motion-free decorrelation value. The regression-based bulk motion subtraction algorithm improved image quality, vessel density measurement repeatability and accuracy of bulk motion noise removal compared to an earlier median subtraction algorithm. Pre-compensation by OCT inter-B-scan registration before OCTA could improve the efficacy in removing visible artifacts at the cost of increased processing time and reduced capillary-background contrast. The regression-based algorithm takes into account the dependence of decorrelation on reflectance signal and hence, it shows potential to reduce the effect of signal strength on vessel identification. Furthermore, it can be used to identify flow voxels without the need for an avascular reference area. Finally, since it corrects flow signal on voxels of the three-dimensional OCTA data set, it might improve capillary flow index quantification accuracy and image quality on projection-resolved OCTA.