The Impact of Image Processing Algorithms on Optical Coherence Tomography Angiography Metrics and Study Conclusions in Diabetic Retinopathy

Abstract

Purpose: The purpose of this study was to evaluate the impact of image processing on quantitative metrics in optical coherence tomography angiography (OCTA) images and study conclusions in patients with diabetes.

Methods: This was a single center, retrospective cross-sectional study. OCTA imaging with the Cirrus HD-OCT 5000 AngioPlex of patients with diabetes was performed. The 8 × 8 mm superficial slab images underwent 4 different preprocessing methods (none, background subtraction [BGS], foveal avascular zone brightness adjustment, and contrast limited adaptive histogram equalization [CLAHE]) followed by 4 different binarization algorithms (global Huang, global Otsu, local Niblack, and local Phansalkar) in ImageJ. Vessel density (VD), skeletonized VD (SVD), and fractal dimension (FD) were calculated. Mixed-effect multivariate linear regressions were performed.

Results: Two hundred eleven scans from 104 patients were included. Of these scans, 67 (31.8%) had no diabetic retinopathy (DR), 99 (46.9%) had nonproliferative DR (NPDR), and 45 (21.3%) had proliferative DR (PDR). Forty-eight of 211 (22.7%) scans had diabetic macular edema (DME). The image processing method used significantly impacted values of VD, SVD, and FD (all P -values < 0.001). On multivariate analysis, the image processing method changed the clinical variables significantly associated with VD, SVD, and FD. However, BGS and CLAHE yielded more consistent significant covariates across multiple binarization algorithms.

Conclusions: The image processing method can impact the conclusions of any given study analyzing quantitative OCTA metrics. Thus, caution is urged in the interpretation of such studies. Background subtraction or CLAHE may play a role in the standardization of image processing.

Translational Relevance: This work proposes strategies to achieve robust and consistent analysis of OCTA imaging, which is especially important for clinical trials.

Introduction

Optical coherence tomography angiography (OCTA) enables noninvasive, high-resolution, and depth-resolved visualization of the retinal microvasculature, and has become more widely utilized in both research and patient care over the past several years. Much research has focused on quantitative metrics derived from OCTA images, such as vessel density (VD), skeletonized vessel density (SVD), and fractal dimension (FD), among others. The values of these metrics have been shown to be altered in retinal vascular disease, such as diabetic retinopathy (DR).^1–5 

The calculation of these quantitative metrics requires image processing of OCTA images which involves binarization, a process that converts the original grayscale image into a black and white image based on a numerical threshold. Many binarization algorithms exist - some are “global,” in which one numerical threshold is determined for the entire image, and some are “local,” where different thresholds are calculated for different areas of the image depending on characteristics of nearby surroundings. There are many global and local binarization algorithms that have been used in various OCTA studies.^2,3,6–10 However, there is no standardization or consensus on how binarization should be done. Furthermore, there are also various preprocessing methods previously used in the literature to adjust brightness and contrast prior to performing binarization.^6,11 Unfortunately, there is no standardization or consensus for the most appropriate way to perform image processing on OCTA images. The lack of standardization prevents comparison of quantitative metrics between studies and has important implications for clinical trials utilizing OCTA. 

Several OCTA platforms, including the RTVue-XR Avanti (Optovue Inc., Fremont, CA), Cirrus HD-OCT 5000 Angioplex, and PLEX Elite 9000 (Carl Zeiss Meditec, Dublin, CA) have recently acquired built-in software to calculate VD. Furthermore, the underlying algorithms or processes used to perform these calculations have not been made publicly available, which prevents comparative assessment of these algorithms as there is currently no ground truth by which to validate the calculation of these metrics. As a result, calculated OCTA quantitative metrics cannot be generalized beyond a single machine. 

Several recent studies have highlighted significant variations in quantitative metrics calculated from OCTA images depending on the image processing procedures used. In a recent single-center study examining 21 eyes from 11 individuals who underwent swept-source (SS-)OCTA processed with different binarization and brightness/contrast adjustment algorithms, Mehta et al. found statistically significant differences between OCTA quantitative measurements from different binarization thresholding methods.¹² In another cross-sectional study of 44 eyes from healthy subjects who underwent OCTA imaging that was processed with 7 well-established algorithms from the literature, Rabiolo, et al. found that different algorithms had poor agreement and poor reliability.¹³ These and other studies raise the question of whether image processing methodology can impact not only the values of quantitative metrics but also the conclusions of any given study. 

Herein, we present a retrospective analysis of 211 OCTA scans from 104 patients with diabetes where we investigated several combinations of preprocessing methods and image binarization algorithms to determine whether and how image processing impacts OCTA-based calculations of VD, SVD, and FD. We also assessed whether the image processing method used could change which clinical variables were statistically significant in multivariate regression models, which would impact the conclusions of the study. 

Methods

Study Design

This was a retrospective, observational study performed at the Yale Eye Center, Yale School of Medicine, New Haven, CT, United States. This study upheld the tenets of the Declaration of Helsinki and was approved by the Institutional Review Board (IRB) at Yale. 

Study Protocol

Patients with a diagnosis of diabetes mellitus from the practices of two retina specialists at the Yale Eye Center (authors K.N. and R.A.) who underwent OCTA imaging from January 1, 2017, to April 1, 2019, were included. The study excluded patients under the age of 18 years, eyes with a history of ocular trauma, and other co-existing chorioretinal disease, such as retinal vascular occlusion, high myopia, and macular telangiectasia. Complete ophthalmologic examinations had been performed at each visit. Clinical and demographic variables including age, sex, stage of DR - no DR (no DR), nonproliferative diabetic retinopathy (NPDR), and proliferative diabetic retinopathy (PDR), last hemoglobin A1c, best-corrected visual acuity, and presence of diabetic macular edema (DME) were recorded. DME was defined as the presence of clinically significant macular edema (CSME) as designated by the Early Treatment Diabetic Retinopathy Study (ETDRS). 

Patients were imaged using the AngioPlex OCTA platform on the spectral-domain Cirrus HD-OCT 5000 (Carl Zeiss Meditec Inc, Dublin, CA). The Cirrus HD-OCT 5000 system scans at 68,000 A-scans per second and has axial and transverse resolutions of 5 µm and 15 µm, respectively. The 8 × 8 mm scan centered on the fovea was utilized because it was the largest field of view available and thus most likely to detect areas of nonperfusion further from the fovea in patients with more advanced stages of DR. Segmentation was automated by the built-in OCTA software version 11.5.3.61246, which uses a multilayer segmentation algorithm. The decision was made to analyze only the superficial capillary network (segmented between the internal limiting membrane and the inner plexiform layer) for the purposes of this study because of the potential for projection artifacts in images from the deep capillary network. Images of poor quality, from low signal strength, media opacity, or motion artifacts, were excluded from image processing. Segmentation was determined to be accurate for all included images. 

Image Analysis

The 8 × 8 mm superficial slab images were automatically exported using the built-in software version 11.5.2.54532 as BMP images and imported to ImageJ (National Institutes of Health [NIH], Bethesda, MD) with dimensions of 1024 by 1024 pixels. Figure 1 summarizes the image processing protocol. All images were converted to 8-bit format. Prior to binarization, images were preprocessed using one of four methods: no preprocessing, background subtraction, brightness adjustment based on grayscale values within the foveal avascular zone (FAZ), and contrast limited adaptive histogram equalization (CLAHE). Background subtraction is performed with the “Subtract Background” command in ImageJ with the default rolling ball radius of 50 pixels. This enables the removal of uneven backgrounds, which, for example, may occur if vitreous floaters or other media opacity result in local reduction of the signal. The technique of brightness adjustment based on the FAZ has been previously reported.² Briefly, the average grayscale value of 3 circles of radius 10 pixels centered around the FAZ is calculated, and all pixel grayscale values less than or equal to this value is set to zero. The default values for CLAHE in ImageJ were used (block size 127, histogram bins 256, maximum slope 3, and with fast analysis). CLAHE is an algorithm that locally enhances contrast but limits amplification of noise, which has also been previously used.^12,14 

Figure 1.

 

Flowchart detailing image processing pipeline from preprocessing steps to binarization. FAZ, foveal avascular zone; CLAHE, contrast limited adaptive histogram equalization.

Images were binarized using one of four algorithms: global Huang, global Otsu, local Niblack, and local Phansalkar in ImageJ. Global binarization algorithms calculate a single grayscale threshold for the entire image. The Huang algorithm calculates a global threshold by minimizing the measure of “fuzziness” using Shannon's entropy function.¹⁵ The Otsu algorithm calculates a global threshold which minimizes the intraclass intensity variance within the foreground and background.¹⁶ Local binarization algorithms calculate different grayscale thresholds for each pixel within the image depending on the image characteristics surrounding the pixel within a certain radius. The “radius” setting on the local binarization algorithms was adjusted to 100 from the default of 15 to eliminate signal from within the FAZ after binarization. The Niblack local thresholding algorithm calculates thresholds based on each pixel's deviation from the mean.¹⁷ The Phansalkar local thresholding algorithm is a modification of the Niblack algorithm to deal with low contrast images.¹⁸ Figure 2 shows a representative image in a patient without DR and resulting images after the above preprocessing and binarization algorithms. Figure 3 shows similar representative images in a patient with moderate NPDR. The VD was calculated as the number of bright pixels divided by the total number of pixels. The FD was calculated with the box counting method using the FracLac plugin in ImageJ. Subsequently, skeletonization of the binarized image was performed, and the SVD was calculated as the number of bright pixels divided by the total number of pixels. 

Figure 2.

 

Representative examples of optical coherence tomography angiography images in a patient without diabetic retinopathy after first applying one of four preprocessing methods, and subsequently one of four binarization algorithms. FAZ, foveal avascular zone; CLAHE, contrast limited adaptive histogram equalization.

Figure 3.

 

Representative examples of optical coherence tomography angiography images in a patient with moderate nonproliferative diabetic retinopathy and capillary nonperfusion after first applying one of four preprocessing methods, and subsequently one of four binarization algorithms. FAZ, foveal avascular zone; CLAHE, contrast limited adaptive histogram equalization.

Statistical Analysis

Univariate mixed-effect linear regression analysis was performed to assess the impact of preprocessing method and binarization algorithm on VD, SVD, and FD to account for the inclusion of more than one eye per patient. Multivariate mixed-effect linear regression analysis controlling for signal strength, age, gender, eye laterality, stage of DR, and presence of DME was performed on VD, SVD, and FD. All statistical analysis was performed using Stata 14 (StataCorp, College Station TX). All graphical analysis was performed using GraphPad Prism 8 (GraphPad, San Diego, CA) and figures were generated using GraphPad Prism 8 and Adobe Illustrator (Adobe Inc., Mountain View, CA). The P values < 0.05 were considered statistically significant. 

Results

Image Characteristics

Three hundred one OCTA images from 124 patients were initially included. Ninety of these scans were excluded for poor signal strength or significant image artifact, resulting in a total of 211 scans from 104 patients analyzed. The mean (standard deviation [SD], range) age of patients was 56.6 years (SD = 14.2, range = 22-88 years), and 61 (58.7%) were women. Of these scans, 67 (31.8%) scans had no DR, 99 (46.9%) had NPDR, and 45 (21.3%) had PDR. Forty-eight (22.7%) scans had DME. More detailed clinical and demographic characteristics are summarized in Table 1. 

Table 1.

 

Demographic, Clinical, and Imaging Characteristics of Optical Coherence Tomography Angiography Scans

 

Quantitative OCTA Metrics

The values of all angiographic metrics (VD, SVD, and FD) were all significantly associated with the choice of the preprocessing method and the binarization algorithm (all P values < 10⁻⁴) as shown in Table 2. Preprocessing with the FAZ adjustment method yielded higher estimates of VD (β = 0.447, standard error [SE] = 0.044) and SVD (β = 0.097, SE = 0.010) when combined with the Huang algorithm, but yielded the lowest estimates of VD, SVD, and FD when combined with the Otsu, Niblack, or Phansalkar algorithms. Conversely, preprocessing with the BGS algorithm yielded the highest estimates of VD when combined with the Huang algorithm alone. For all other binarization algorithms, the CLAHE preprocessing method yielded the highest estimates of VD, SVD, and FD (see Table 2). 

Table 2.

 

Values of Quantitative Angiographic Metrics by the Preprocessing Method and the Binarization Algorithm

 

Examining the angiographic metrics by stage of DR, there was a trend of decreasing VD, SVD, and FD with increasing severity of DR. This trend was present regardless of the combination of preprocessing and binarization algorithm used (Fig. 4). 

Figure 4.

 

Violin plots of vessel density (VD), skeletonized vessel density (SVD), and fractal dimension (FD) across the dataset based on preprocessing method and binarization algorithms used. DR, diabetic retinopathy; NPDR, nonproliferative diabetic retinopathy; PDR, proliferative diabetic retinopathy; BGS, background subtraction; FAZ Adj, foveal avascular zone adjustment; CLAHE, contrast limited adaptive histogram equalization.

There was a positive association between signal strength and all OCTA metrics that was statistically significant (P values < 0.001), so signal strength was adjusted for in all subsequent multivariate analyses. 

Multivariate Regression Analyses for Stage of DR and DME

The variables that were significant on multivariate mixed-effects linear regression models varied based on the selection of preprocessing and binarization algorithms (P values in Table 3). The values of each coefficient and 95% confidence intervals are shown in Supplementary Table S1. For instance, age was a significant predictor of decreasing VD for all combinations of preprocessing and binarization except when using the Huang binarization algorithm without preprocessing or with the FAZ adjustment method. Similarly, presence of PDR significantly predicted decreased SVD and FD for nearly all combinations of preprocessing and binarization, except when using the Huang and Otsu binarization algorithms without preprocessing or FAZ preprocessing. The presence of NPDR was consistently not a significant predictor of VD or SVD across all image processing algorithm choices except when using the Phansalkar binarization algorithm with CLAHE. The significance of DME status varied depending on the combination of preprocessing and binarization algorithms. When using the Huang binarization algorithm with the FAZ adjustment method, no clinical variables were statistically significant. When comparing all preprocessing methods, the use of background subtraction most consistently resulted in the same clinical variables being statistically significant regardless of which binarization algorithm was subsequently used. In addition, CLAHE also yielded mostly consistent results in this respect. 

Table 3.

 

Table and Heat Map of P Values From Multivariate Mixed-Effects Linear Regression Models Across Preprocessing and Binarization Algorithms for Age, Stage of Diabetic Retinopathy, and the Presence of Diabetic Macular Edema

 

Discussion

In this retrospective analysis of 211 OCTA scans from 104 patients with diabetes, we found that both the choice of the preprocessing method and the image binarization algorithms significantly impacted the values of quantitative angiographic metrics of VD, SVD, and FD. Furthermore, the choice of the preprocessing method and the binarization algorithm also impacted which clinical variables were significantly associated with these metrics. Thus, the way image processing is done can potentially alter the conclusions of any given study. However, we also found that if background subtraction was done prior to binarization, the clinical variables that were found to be statistically significantly associated with the quantitative metrics were mostly consistent regardless of the binarization algorithm used. Fairly consistent results were also obtained using CLAHE prior to binarization. 

There is a wealth of OCTA literature in DR, and numerous studies that have reported reduced VD, SVD, and FD in patients with DR.^2–4,10 It is also known that older age is associated with decreased VD and SVD.^5,19 This is consistent with our findings that PDR and older age were most frequently significantly associated with decreased values of VD, SVD, and FD, which served as somewhat of a positive control in our study to evaluate which combinations of preprocessing and binarization algorithms might be most accurate. 

Other studies have examined in more detail the issue of image processing in OCTA and its effects on measurements of quantitative metrics. For instance, in a study of 15 healthy control subjects and 30 patients with DR using 3 × 3 mm OCTA images, Borrelli et al. showed that using global binarization algorithms yielded significantly lower measurements of quantitative OCTA vascular parameters compared to local binarization algorithms.²⁰ They found this to be particularly true in eyes with DME, which they hypothesized was due to thresholding of “pseudoflow” in the macula by local binarization algorithms. This was not necessarily the case in our study, as whereas the local binarization algorithm Phansalkar did yield higher measurements than the other binarization algorithms, the Niblack local algorithm did not yield consistently higher measurements than the Huang global algorithm. In addition, DME was not consistently associated with higher quantitative measurements in local binarization algorithms. The most likely explanation for this discrepancy is the selection of the radius when using local binarization algorithms. In the study by Borrelli et al. and in other previous studies, the default radius of 15 pixels was chosen, whereas, in our study, we chose 100 pixels.^14,20 The reason is that when the default value of 15 pixels is used, the signal within the FAZ is inappropriately thresholded. We considered thresholded signal within the FAZ inaccurate and thus increased the radius until the thresholded signal no longer appeared in the FAZ, which occurred for most images at around 100 pixels (approximately 780 µm), which is why this value was chosen. The smaller the radius, the more the threshold will vary by location within the image. Conversely, the larger the radius, the more the local algorithm resembles a global algorithm as there will be less threshold variability. In general, more signal will be thresholded with a smaller radius, thus increasing the values of the quantitative metrics calculated.²¹ This underscores the importance of all details and parameters when it comes to image processing methodology. Further work is needed to determine the optimal radius for local binarization algorithms which will likely depend on the size of the OCTA scan, the depth of the scan, as well as the resolution of the image itself in pixels. We speculate that the optimal value will likely be on the order of the diameter of the FAZ, as in this study. 

In a cross-sectional study by Rabiolo et al. using 6 × 6 mm OCTA images from the PLEX Elite 9000 (Zeiss), 7 previously published binarization algorithms were investigated, including the proprietary algorithm from the Zeiss ARI (Advanced Retinal Imaging) Network, and the authors found that different algorithms were not directly interchangeable because they produced significantly different quantitative values.¹³ These authors included healthy eyes, eyes with DR, and eyes with glaucoma, and found that eyes with DR and glaucoma had lower vessel densities regardless of binarization algorithm. They, however, did not study the influence of preprocessing methods prior to binarization or local binarization algorithms. 

A smaller single-center study by Mehta et al., utilized 3 × 3 mm OCTA images of healthy eyes from the PLEX Elite 9000 (Zeiss) and processed the images with different preprocessing methods as well as global and local binarization algorithms. They also found statistically significant differences between quantitative measurements from different image processing methods, and that brightness or contrast adjustments interacted with binarization algorithm to produce different results.¹² Furthermore, the same group found poor repeatability for a wide range of preprocessing methods and binarization algorithms, further urging caution when analyzing quantitative OCTA parameters.¹⁴ 

In another study utilizing the PLEX Elite 9000 by Hong et al., 3 × 3 mm scans were compared to 12 × 12 mm scans in healthy patients, and images were processed using different vessel enhancement filters before using a single binarization method to calculate perfusion density and VD.²² They found that the repeatability of these quantitative metrics varied depending on the vessel enhancement filter used, and only moderate correlation at best between the 3 × 3 mm and 12 × 12 mm scans. These authors also applied similar methods to even larger 15 × 9 mm scans and found that the vessel enhancement filter used affected the repeatability of the VD measurements, further emphasizing the importance of being consistent in image processing.²³ These studies, however, were performed on nondiabetic patients and did not explore the effect of using different binarization algorithms. 

One advantage to using local binarization algorithms over global binarization algorithms is the potential to threshold the signal regardless of regional variations of brightness or intensity within the image. However, as described above, local binarization algorithms may also have the tendency to threshold the signal that is “noise” and not truly representative of flow signal. This is seen in thresholding of the signal within the FAZ, which was mitigated in this study by increasing the radius parameter. Regardless, in the absence of an accepted ground truth for the retinal microvasculature, it is difficult to determine whether one method is in fact superior. 

With regard to preprocessing methods, the goal is to modify brightness and contrast in such a way that enables more accurate binarization. The FAZ adjustment preprocessing method was found to yield lower estimates of VD, SVD, and FD than most other algorithm combinations, whereas preprocessing with CLAHE yielded the highest estimates for VD, SVD, and FD across most algorithm combinations, especially when paired with the Phansalkar algorithm. Notably, CLAHE also appeared to decrease the variability of the quantitative metrics across binarization algorithms, which is consistent with the previous study by Mehta et al. CLAHE also yielded relatively consistent results in terms of which clinical variables were found to be significantly associated with decreased angiographic metrics.¹² Despite these results, we noticed that CLAHE does appear to alter the image histogram more than other preprocessing methods, and thus caution should be used when considering this method. 

We found that the preprocessing method of background subtraction yielded the most consistent results in terms of which clinical variables were significantly associated with quantitative metrics regardless of the binarization algorithm used. The background subtraction command in ImageJ uses a “rolling ball” algorithm to calculate a local background value for every pixel within a specified radius, and then subtracts that value from the original image.²⁴ It is advised that the specified radius be larger than the largest object that is not part of the background. In the context of superficial OCTA images, this would be the width of the largest retinal vessel in the image. In this study, the default radius of 50 pixels was used but it is unclear whether there is an optimal value. Background subtraction, in effect, can compensate for variable background brightness throughout the image that may be caused by media opacity and shadowing, producing a more uniformly bright image that subsequently could potentially be more accurately binarized. 

This study has several limitations which may influence its interpretation and generalizability. First, as with all retrospective studies, the present work is limited to correlative, rather than causal, conclusions about the associations between DR and vascular parameters. Second, although we explore a range of possibilities for preprocessing and binarization of OCTA images, our analysis is not exhaustive, and caution should be used when applying these results to other algorithms not described in this work. Furthermore, the current study only focuses on 8 × 8 mm superficial vascular plexus OCTA images of patients with diabetes. Although we did not perform image averaging to investigate repeatability and reproducibility of image acquisition and perhaps improve image quality, we did control for signal strength as a measure of image quality in our analyses. In addition, the Zeiss Cirrus machines at our institution are not equipped with the AngioMetrix quantitative analysis software as it is not approved for use in the United States by the US Food and Drug Administration (FDA), and thus we were unable to perform a quantitative comparison with angiometric values provided by the device itself. We also did not investigate other quantitative OCTA metrics, such as the FAZ area because these are less dependent on image processing. Despite these limitations, this work provides clear evidence of the impact of image processing methodology on vascular metrics derived from OCTA images and in the clinical variables that are found to be significantly associated with these vascular metrics. Further work may aim to automate and standardize the processing and analysis of OCTA images. 

In conclusion, we demonstrate the impact of OCTA image analysis methodology on any given study's results and urge caution when interpreting such studies. There is a need for standardization of image analysis methodology to ensure accurate and consistent results with implications for clinical trials utilizing OCTA for monitoring of disease progression using quantitative OCTA metrics. We also believe there is a need for detailed documentation of image processing methodology across all studies and transparency in proprietary software so that these algorithms can be faithfully evaluated against each other, and especially against an accepted ground truth if one emerges in the future. In the meantime, we suggest considering the use of background subtraction or potentially CLAHE as a preprocessing step, as in our study, this led to more consistent results robust to the binarization algorithm used. We also suggest that more than one binarization algorithm be used in any given study to ensure that the results remain consistent and robust across multiple image analysis methods. Thoughtful discussion is needed on what should be considered the most appropriate way to perform image processing on OCTA images. 

Acknowledgments

This work has been previously presented at the Retina Society Annual Meeting, 2021 in Chicago, IL, and the American Society of Retina Specialists Annual Meeting, 2021 in San Antonio, TX. 

Disclosure: I.G. Freedman, None; E. Li, None; L. Hui, None; R.A. Adelman, None; K. Nwanyanwu, None; J.C. Wang, None