Fully automated geometric feature analysis in optical coherence tomography angiography for objective classification of diabetic retinopathy

: This study is to establish quantitative features of vascular geometry in optical coherence tomography angiography (OCTA) and validate them for the objective classification of diabetic retinopathy (DR). Six geometric features, including total vessel branching angle (VBA: θ ), child branching angles (CBAs: α 1 and α 2), vessel branching coefficient (VBC), and children-to-parent vessel width ratios (VWR1 and VWR2), were automatically derived from each vessel branch in OCTA. Comparative analysis of heathy control, diabetes with no DR (NoDR), and non-proliferative DR (NPDR) was conducted. Our study reveals four quantitative OCTA features to produce robust DR detection and staging classification: (ANOVA, P<0.05), VBA, CBA1, VBC, and VWR1.

measuring th ophthalmolog used to predic major limitat imaging for q qualitative inf limited ability retina [21]. T geometric dis the geometric the geometric progression. A (OCTA) whic resolution eve

Materials
The detailed illustrated in branch triang algorithm use Fig tification. ging using a fully automated algorithm. These features can be divided into two subcategories, vessel angle-based and vessel width-based geometric features. Vessel angle-based features measure the degree of bifurcation and relative angle change among the parent and child vessel branches. The first feature is vessel branching angle (VBA: θ), which measures the overall degree of bifurcation between the child branches. The second and third features are the child branching angle 1 (CBA1: α1) and child branching angle 2 (CBA2: α2), which measures the direction of branching of each child vessel. The child angles can provide information relating to the deviation of direction of child vessels with respect to the original parent vessel direction.
Vessel width-based features measure the structural change of the vessels as a result of branching. The first vessel-width feature is vessel branching coefficient (VBC), which measures the effect of branching on vessel area between the parent, first (larger) child, and second (smaller) child. Similarly, the second and third vessel-width features are vessel width ratio 1 (VWR1) and width ratio 2 (VWR2), which examines the change in width between individual child and parent vessels.

Branching triangle identification
In this step, the processed vessel map is used to extract the vessel skeleton, branch points, and end points. Figure 2(B) illustrates the steps for locating branch points and corresponding endpoints. First, the vessel skeleton map is extracted from the processed vessel map. Next, the vessel branchpoints are determined from the vessel skeleton using morphological functions ( Fig.  2(B1)). The branchpoints indicate where vessel branching occurs. Furthermore, we can utilize the branchpoints to create a mask to extract the region of interest ( Fig. 2(B2)). In this study, the region of interest is referred to as branch triangles. The method of masking starts by dilating each branchpoint. Based on empirical trials, a dilation radius of 5 pixels is large enough to encompass the region of interest, and small enough to not overlap with adjacent branchpoints. The newly created mask is then multiplied with the skeleton map to generate an image containing the individual branch triangles. The branch triangle visualization is referenced as the colored pixels in Fig. 2(B2). In the last step the endpoints of each branch triangles are extracted into a separate image (Fig. 2(B3)). The preprocessing and branch triangle identification steps results in four images, i.e., the processed vessels, vessel skeleton, branchpoints, and endpoints, used for following feature measurement.

Quantification of features
The inputs into the fully automated algorithm are the processed vessel, vessel skeleton, branchpoint, and endpoint images. The algorithm is illustrated in Fig. 1(b) and runs iteratively for each branchpoint in the image. The algorithm is designed to quantify each branchpoint with exactly 3 endpoints. For the point that is crossed by two large vessels, the corresponding number of endpoints will be 4, and therefore is excluded from quantitative branchpoint analysis.
For a branchpoint, the first step is the branch triangle endpoint identification by a window method illustrated in Fig. 2(B3). For this process, the window area is centered at the branchpoint, and the window is a square whose width is 15 pixels. The width size is dependent on the dilation radius used in the branch triangle identification step. If the window is too small, then none of the endpoints would be detected. Similarly, if the width is too large, then the endpoints of other branches would be detected. The window method records the coordinates of each endpoint of the branch triangles.
Next the recorded coordinates of the endpoints and branchpoints are used to create vectors of the vessels, resulting in a set of parent and two child vectors. From these three vectors, three angles, A and B (Fig. 2(C1)), and θ ( Fig. 2(C2)) are quantified ( Fig. 1(b)). Based on the physiology of vessel branching, the smallest one of these three angles is the VBA (θ). Using geometric identities, CBA1 (α1) and CBA2 (α2) are calculated from angle B and A, respectively, of the vessel, where CBA1 is larger than CBA2. We can then use the calculated angles (θ, α1, α2) to identify the parent and child vessels. The two vectors/endpoints that form the branching angle are identified as the child vessels, and the remaining vector/endpoint identifies the parent vessel.
By identifying the parent and child vessels, the width of each vessel is measured by using another window method. For this process, the window area is centered around each endpoint (Fig. 2(D1)). The coordinates of the window area are used to take the ratio of the sum of the of pixels in the processed vessel map (Fig. 2(D2)) to sum of the pixels in the vessel skeleton map (Fig. 2(D3)). This ratio defines the width of the vessel segment. The window method is repeated for all parent and child vessels.
The values of the widths are then used to calculate the width parameters, VBC, VWR1, and VWR2. The width parameters are summarized in equations Table 1

Statistical analysis
The statistical analyses were performed using MATLAB (MathWorks, Natick, MA) and Microsoft Excel (Microsoft Corporation, Redmond, WA). All quantitative geometric features for different groups were compared using one-way, multi-label analysis of variance (ANOVA) to compare difference in mean values of the features among multiple groups. Following ANOVA, one versus one comparisons of the features between control, NoDR and NPDR stages (mild, moderate and severe) were performed using a Student t-test, with statistical significance was defined with a P value < 0.05.

Results
Our OCTA image database consisted of 120 OCTA images from 60 NPDR patients, 24 OCTA images of 17 NoDR patients, and 40 images from 20 control subjects. We excluded 6 eyes due to poor OCTA images and another 4 eyes due to presence of other retinal diseases. Among the NPDR OCTAs, there were 20 sets each of mild, moderate, and severe NPDR (40 OCTAs for each stage). There was no statistical significance in the distribution of age, sex and hypertension between control and DR groups. (ANOVA, P = 0.14; chi-square test, P = 0.11 and P = 0.32, respectively). The demographics of our studied subjects are detailed in Table 2.
The individual t-test comparison for VBA (θ) is shown in Table 4. In the case of VBA (θ), control groups can be significantly differentiated from NoDR, mild, and moderate groups using individual t-test (Student's t-test: P<0.05 for control vs NoDR, mild, and moderate NPDR). NoDR groups can be highly significantly differentiated from mild, and moderate groups (Student's t-test: P<0.001 for NoDR vs mild, and moderate NPDR), and mild can be significantly differentiated from severe groups (Student's t-test: P<0.05 for mild vs severe NPDR). The individual t-test comparison for CBA1 (α1) is shown in Table 5. For CBA1 (α1), NoDR groups can be significantly differentiated from mild, and moderate groups (Student's ttest: P<0.01 for NoDR vs mild, and moderate NPDR), while severe groups can be significantly differentiated from mild, and moderate groups (Student's t-test: P<0.05 for severe vs mild, and moderate NPDR). The individual t-test comparison for VBC is shown in Table 6. In case of VBC, control groups can be significantly differentiated from NoDR and severe groups (Student's t-test: P<0.01 for control vs NoDR, and severe NPDR). No other groups can be differentiated from each other. The individual t-test comparison for VWR1 is shown in Table 7. For VWR1, control groups can be significantly differentiated from NoDR, mild, and severe groups (Student's ttest P<0.05 for control vs NoDR, and mild NPDR, and P<0.01 for control vs severe NPDR). No other groups can be differentiated from each other.

Discussion
In summary, we developed a fully automated method to quantify geometric features, including VBA, CBA1, CBA2, VBC, VWR1, and VWR2. Comparative OCTA analysis of healthy control, diabetes with NoDR and NPDR subjects were conducted.
As stated by Murray's principle, optimal blood flow transport is achieved with minimum energy and maximum vascular diffusion into surrounding tissue [3]. Changes to the vascular structures caused by DR and other diseases would decrease optimal transport of blood and nutrients. Therefore, quantitative geometric analysis can provide biomarkers for objective DR detection, progression monitoring, and treatment assessment. In this study, we quantified two categories of geometry, i.e., the angle (VBA, CBA1, CBA2) and width-based (VBC, VWR1, and VWR2) features, in OCTA.
In principle, the angle-based features can quantify abnormal changes in blood transport optimality in DR. Using traditional fundus photography, previous studies has observed widening of VBA associated with DR [15,25,26]. However, fundus photography has limited resolution and contrast to disclose subtle changes, particularly in the foveal area where the capillary structures are small, at early onset of mild DR [27]. Moreover, previous studies used manual or semi-manual marking software to analyze the geometric features, which was time consuming with labor intensive involvement of experienced ophthalmologists [15,25,28]. By providing depth-resolved imaging capability, OCTA allows high resolution visualization of micro vasculatures in the retina [29,30]. Quantitative features, such as tortuosity, vessel density, etc., in OCTA have been used for evaluating DR [31,32], SCR [33,34], etc., However, quantitative analysis of vessel angle analysis is not yet established. In this study, we demonstrated quantitative analysis of VBA automatically. To further quantify the bifurcation, the automated algorithm could also quantify individual angles of the child vessels, i.e., CBA1 and CBA2, to evaluate the orientation changes from the parent vessel direction. We observed that VBA (θ), CBA1 (α1), and CBA2 (α2) increased. This indicates that both vessels contribute to the widening of the total VBA.
Similarly, width-based features, such as VBC, also reflect transport optimality in retinal vasculature. Vessel width is related to the vessel cross sectional area and volume, and changes to the vessel width will affect the flow of blood and nutrients. Quantitative OCTA analysis of vessel caliber has been demonstrated [31,32]. However, comprehensive analysis of vessel branching properties is not yet validated. In this study, we verified VBC as a sensitive biomarker to DR detection and staging classification. By using ratio parameter VBC, the effect of image size and scale on quantitative measurement can be minimized. We also measured VWR1 and VWR2, width ratios of the individual child vessels to the width of the parent vessel, to provide further information such as contraction or dilation of child vessels. Based on the report by Adam et. al. [35], symmetrical bifurcations should have a branching coefficient of 1.26, and in case of non-symmetrical bifurcations, the VBC decrease from this optimal value. In our study, we observed that healthy control subjects had an average VBC of 1.20, close to the theoretical optimal value. We also observed that the VBC decreases in DR patients compared to healthy control patients. Similarly, we observed decreased VWR1 and VWR2 values in DR patients. We speculate the decreased vessel widths and increased vessel branching angles are possibly due to remodeling related to ischemia and neovascularization. This remodeling could increase the shear stress at bifurcations resulting in increased vessel branching angles. Further functional imaging study of retinal hemodynamics may provide insights to a better understanding of these geometric parameter changes.
It was observed that for NPDR, VBA widens from mild to moderate groups, but then narrows from moderate to severe (Table 3). This observation is consistent to previous study with traditional fundus photography [26]. Interestingly, we observed that VBA decreases from control to NoDR (Table 3), then increases from NoDR to moderate, and then decreases once more from moderate to severe. This phenomenon might reflect an auto-regulatory response due to increased metabolic demand in NoDR groups. Jiang et. al. [36] recently observed a similar trend in retinal thickness in NoDR OCT. These unique OCTA features might reflect subtle vascular abnormalities that cannot be revealed in traditional fundus photography.

Conclusion
In conclusion, an automated algorithm was developed to quantify six geometric features, including VBA (θ), CBA1 (α1), CBA2 (α2), VBC, VWR1, and VWR2. Comparative analysis of healthy control, diabetes with NoDR and DR subjects revealed that the VBA, CBA1 (α1), VBC, and VWR1 showed significant difference among control, NoDR and NPDR groups. The study demonstrates the potential of using quantitative geometric features as objective biomarker for DR detection, progression monitoring, and treatment assessment.