2D alpha-shapes to quantify retinal microvasculature morphology and their application to proliferative diabetic retinopathy characterisation in fundus photographs

The use of 2D alpha-shapes (α-shapes) to quantify morphological features of the retinal microvasculature could lead to imaging biomarkers for proliferative diabetic retinopathy (PDR). We tested our approach using the MESSIDOR dataset that consists of colour fundus photographs from 547 healthy individuals, 149 with mild diabetic retinopathy (DR), 239 with moderate DR, 199 pre-PDR and 53 PDR. The skeleton (centrelines) of the automatically segmented retinal vasculature was represented as an α-shape and the proposed parameters, complexity (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Op\alpha }_{min}$$\end{document}Opαmin), spread (OpA), global shape (VS) and presence of abnormal angiogenesis (Gradα) were computed. In cross-sectional analysis, individuals with PDR had a lower \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${Op\alpha }_{min}$$\end{document}Opαmin, OpA and Gradα indicating a vasculature that is more complex, less spread (i.e. dense) and the presence of numerous small vessels. The results show that α-shape parameters characterise vascular abnormalities predictive of PDR (AUC 0.73; 95% CI [0.73 0.74]) and have therefore potential to reveal changes in retinal microvascular morphology.

. Such variation in previously reported FD have been attributed to lack of standardisation in FD calculation 17,18 . FD and it's calculation can be affected by: (1) the accuracy of vessel segmentation; (2) skeletonization, locating the centreline of each vessel; (3) camera resolution and image quality, which directly affects (1) and (2) 17 ; (4) the FD algorithm and the choice of its parameters, such as the random sampling of points on the vessel map; and (5) the computational cost 19 , which may be excessive in a clinical setting.
In this study, we explore an alternative method of quantifying the global morphology of the retinal vasculature, using α-shapes. The α-shape is a well described geometric concept used to represent the boundary of a shape in 2D or 3D 20 . It has been applied previously to molecular shape characterisation 21 , brain tumour measurement in magnetic resonance imaging (MRI) 22 , cosmic web topology 23 , surface reconstruction 24 and predicting species distributions 25 . Unlike FD, α-shapes are independent from sampling space, do not make assumptions of self-similarity, are computationally inexpensive and unaffected by parameter choice as no random sampling is conducted and there is only one α-shape for a point set 20 . We applied α-shapes to the task of characterising PDR, with a focus on deriving objective measurements that have a clinically interpretable output.

Methods
Materials and data. MESSIDOR is a publicly available dataset commonly used to evaluate performance of automatic DR detection systems and consists of 1187 fundus images acquired by a 45° FOV Topcon TRC NW6 non-mydriatic camera 26 . Of these 547 are grade 0 (no DR), 149 grade 1 (mild DR), 239 grade 2 (moderate), and 252 grade 3 (severe DR, consisting of pre-proliferative and PDR cases). Images are stored in TIF format and have three sizes: 1140 × 960 pixels, 2240 × 1488 pixels and 2304 × 1536 pixels. Segmentations and PDR image labels from Orlando et al. 13 were accessed at https:// github. com/ ignac iorla ndo/ fundus-fract al-analy sis. In 13 , the authors standardised all images to 909 × 909 pixels to approximately match the calibre of the vessels appearing in images in the Digital Retinal Image for Vessel Extraction (DRIVE) dataset 27 . Images were then segmented using a fully connected conditional random field model 28 . Vessel skeletonization was performed using an iterative thinning algorithm that extracts the centre line of the vessels (1-pixel wide; curved lines), resulting in a vessel skeleton 29 . PDR labels were obtained by relabelling of DR grade 3 by two ophthalmologists whereby 53 images of grade 3 were assigned to the PDR grade (with neovessels) 13 . It should be noted that the relative number of images per DR group differ from 13 due to the corrections outlined on the MESSIDOR project website which involved removal of 13 duplicate images and corrections of DR grades. We calculated FD using a multifractal approach and a generalised sandbox method, reported elsewhere 5 .
For assessing accuracy of α-shape descriptors and FD we used the publicly available high resolution fundus (HRF) dataset (available here: High-Resolution Fundus (HRF) Image Database (fau.de)). The HRF dataset contains 15 images of healthy people and 15 images of people with diabetic retinopathy (DR grades for these images are unavailable) with manual annotations of blood vessels. Images are 3504 × 2336 stored in JPEG format. Images were standardised to 1525 × 1017 pixels to approximately match the calibre of the DRIVE dataset. Images were segmented using the same procedure as the MESSIDOR dataset 13 .
2D α-shapes. First, the retinal vasculature from a colour fundus photograph is skeletonised to obtain a map of centrelines of the blood vessels. Then a Delaunay Triangulation ( DT ) is computed from the points S of the skeleton ( Skel ) (Fig. 1a,b). The α-shape can be considered as a subset of DT , whereby Delaunay simplices-i.e., triangles, edges, vertices-are removed to refine the boundary of the shape. Refinement is achieved by rolling a disc with radius α along each point in S , with an arbitrary starting point. If the next point touches the disc-i.e. the two points lie on the boundary of the disk with radius α and contains no other point in S-this becomes the new pivoting point and the corresponding edge of DT is collected. All the Delaunay simplices that meet this criterion forms the α-shape. For α = 15 (pixels) (Fig. 1c), few points meet this criterion and forms an α-shape that includes a small number of DT edges. As α is increased, more and more edges of DT are collected (Fig. 1d,e). Eventually, there is a minimum value of α that encompasses all points of the vasculature giving a single region and a representation of the vessel shape (Fig. 1f). Finally, very large α equals DT as all edges are eventually included and the shape is not refined. We considered the optimum α-shape of the retinal vasculature (Skel α ) to be the minimum value of α that encompassed all points in the skeleton giving only a single region (Fig. 1f).
Computing 2D α-shape descriptors. For each image in the dataset the α-shape was computed (alphaShape function, Matlab r2019b, Mathworks, USA) applied to the vessel skeleton. We define complexity as the property of a vessel skeleton being composed of many parts (i.e., many branches and vessel paths); a complex vasculature would require more refinement (i.e., smaller α) than a less complex vasculature. Therefore, we extracted the minimum value of α that encompassed all points in the skeleton giving only a single region, termed Opα min , as a descriptor of the complexity of the vascular network. However, used alone Opα min would not completely describe shape complexity, as a vessel skeleton may have a small Opα min due to vessels being close together rather than complex with lots of branches. To model how the vessels are spread around the retina, the area of the optimum α-shape (in pixels), OpA, was computed. Using complexity and spread, we can define vessel shape (VS), We evaluated the sensitivity of VS against changes in skeleton shape by incremental morphological erosion of vessel skeletons. Each skeleton was eroded to approximately 30% of the original size (using the Matlab function www.nature.com/scientificreports/ bwmorph with parameter 'spur') to remove the ends of the vessels (20 iterations) 29 (Fig. 2). A linear trend was observed and quantified by calculating the gradient of the robust regression line that fits the percentage of pixels in each erosion against VS , termed Grad α , to give an additional descriptor of vessel shape (Fig. 3). Intuitively, if the vasculature consists of many small branches, VS will decrease at a lower rate. This is because the main vessels that create the overall shape will be preserved in early iterations of skeleton erosion. Likewise, if the vasculature consists of fewer small branches,VS will decrease at a higher rate as the main vessels are eroded in early iterations. We therefore hypothesised that Grad α could suggest the level of fragmentation of the vasculature, hence the presence of many small vessels typical of neovascularisation.
Experiments and statistical analysis. Three experiments were conducted to address the following aims: (1) to assess the level of significance of the differences between all grades of DR (grades 1 to 3 and PDR); (2) evaluate the contribution of each α-shape descriptor and FD as features for PDR classification; and (3) assess the accuracy of α-shape descriptors and FD values. Statistical analysis was conducted in R (version 3.6.0) and results were considered statistically significant at the significance level 0.05. In experiment (1), all groups (Grade 0 to 3 and PDR) were compared using Kruskall-Wallis test and p-values adjusted for multiple comparisons using Benjamani-Hochberg (data was not normally distributed with homogenous variance; Levene test, p > 0.05). Dunn post-hoc analysis was performed to assess which independent variable differed and p-values corrected for multiple comparisons using Bonferroni.
In experiment (2), α-shape descriptors as features for PDR characterisation were modelled using regularised logistic regression by grouping the data into PDR vs all other DR grades (grades 0 to 3). This was to identify features that could be used to distinguish neovascularisation (a defining feature of PDR) compared to other DR grades. Implementation of regularised logistic regression was performed using the R package glmnet (version 3.6.0) that performs model estimation and feature selection simultaneously 30,31 . We experimented with different www.nature.com/scientificreports/ strengths of regularisation (λ) and a choice of regularisation ( α ∈ [0, 1] ). Lasso regularisation (L 1 regularisation; α = 1), ridge regularisation (L 2 regularisation; α = 0) and a mixture of the two, elastic-net (α = 0.5). Lasso tends to discard correlated variables and select a single feature; ridge tends to equally weight features and shrink coefficients towards each other and elastic-net will exclude groups of highly corelated variables. As the α-shape descriptors are expected to be highly correlated, we applied lasso and elastic-net to perform feature selection. Given that PDR cases consisted of 4.4% of the dataset and different sampling of the training data will give rise to different features being selected, we performed bootstrap analysis. A total of 500 bootstrap trials was performed and for each bootstrap trial the data was divided to keep proportions of classes equal (40 grade 0-3; 40 PDR) and model fitting was performed using fivefold cross-validation. Area Under the Curve (AUC) was used to assess model performance in each fold where the average AUC over all bootstrap trials is reported. fivefold cross validation simultaneously selects a value of λ where the highest AUC on the validation set is achieved (λ min ) and the highest AUC on the validation set within one standard error of λ min (λ 1SE ). Feature importance was assessed by counting the number of times each coefficient reached a non-zero weight in each bootstrap. As λ min and λ 1SE will return different numbers of features (i.e. λ 1SE will contain less features) we assessed feature importance for both values of λ 32 .  Evaluating the sensitivity of VS against changes in skeleton shape. Incremental removal of pixels (expressed as a percentage of the number of pixels of the original skeleton) against the candidate VS parameter for images in Fig. 2. As the skeleton is degraded the VS parameter also decreases, indicating that it is modelling change in vascular properties.  (3), the Dice similarity coefficient (DSC) was used to evaluate performance of the automatic segmentation. To investigate variation in α-shape descriptors and FD measures obtained from the manual and automatic segmentation we use relative error (manual annotation as reference). Relative error is obtained by |AV-MV|/AV, where AV is the value of α-shape descriptors or FD obtained from the manual segmentation and MV is the value of α-shape descriptors or FD obtained from the automatic segmentation. To test whether measures obtained manually and automatically were correlated we used Pearson correlation coefficient test. Healthy and DR α-shape descriptors and FD measures were not normally distributed and therefore compared using Wilcoxon rank sum.

Association of α-shape descriptors and FD with PDR. Images with PDR showed a significantly lower
OpA when compared to Grade 0 (p = 0.001), Grade 1 (p = 0.01), Grade 2 (p = 0.02) and Grade 3 (p = 0.05), indicating a vasculature with less spread. A significantly lower Grad α was observed in PDR when compared to Grade 0 (p = 0.002), Grade 2 (p = 0.02) and Grade 3 (p = 0.02) suggesting that PDR contained numerous small vessels. For VS there was no evidence of significant difference between groups following post-hoc testing and p-value adjustment. PDR exhibited a lower FD when compared to all other groups, indicating a vasculature that is less space filling (p < 0.001) ( Table 1). Table 2 shows the number of times each feature reached a non-zero coefficient and average bootstrap performances. Gradα was the strongest feature for PDR classification. The feature that contributed the least for PDR classification was Opα min . Highest performances was achieved using all α-shape descriptors and FD. Table 3 shows the mean relative errors for α-shape descriptors and FD with the manual segmentation as reference. Higher performances of the automatic segmentation were achieved in healthy individuals (DSC = 0.81870) compared to individuals with DR (DSC = 0.58530) demonstrating a high number of erroneous automatic segmentations in the DR group. Alpha-shape parameters exhibited mean relative errors between 10 and 43%, where a significant correlation was observed in Grad α from Healthy individuals and FD from Healthy and DR. Greatest maximum relative error was observed in OpA measured from an individual with DR (manual OpA = 38; automatic OpA = 90; Max = 100%; DSC = 0.58528) who had large number of microaneurysms, haemorrhages and exudates that obstructed the vasculature and resulted in a sparse and fragmented segmentation. DR exhibited a higher Opα min , lower OpA (manual and automatic) compared to Healthy, indicating a vasculature that was less space filling, but was not significant between groups in measures derived from the manual or automatic segmentation. Grad α was lower in DR compared to healthy individuals indicating the presence of small vessels but was not significant in measures derived from the manual or automatic segmentation. FD exhibited a 0.9% relative error in Healthy and a 1.4% error in DR individuals and was correlated to FD derived from manual segmentations. FD was significantly lower in Healthy (p = 0.01) and DR (p = 0.01), indicating a less space filling vasculature.

Discussion
We implemented α-shapes as a method to quantify retinal microvasculature morphology from fundus imaging and investigated its application to the task of PDR characterisation. A significantly lower FD was observed in PDR compared to all other DR grades suggesting that the vasculature is less space filling in PDR. We did not observe the increase in FD of PDR eyes as reported in 13 , likely due to the removal of duplicates, alteration of DR grades (as specified on the MESSIDOR project website) and parameter choice (e.g. random sampling on the vessel map), highlighting a challenge with the repeatability of FD analysis. The α-shape from which the features proposed here are extracted would not be affected by parameter choice as no random sampling is conducted and there is only one α-shape for a point set 20 .
We experimented with Opα min as a measure of complexity and OpA as a measure of vessel spread and observed that PDR had a significantly lower OpA compared to all other DR grades and a lower Opα min compared to grade 0. This suggests that this the vasculature was dense (i.e. the vessels occupy a smaller area and consists of many 'parts'). OpA was determined to be the most distinguishable feature of PDR and was the most important feature for PDR classification for lasso and elastic-net (when Grad α and FD were not included as features for classification), suggesting that a decreased vessel spread could translate to a predictive marker of progression to PDR.
When Grad α was introduced into the model (without FD) this became the most important feature. PDR had a lower mean Grad α compared to all other grades suggesting the vasculature contains numerous small vessels indicating the presence of neovascularisation. However, the limitation to this parameter is that the presence of numerous small erroneous segmentations of PDR related pathology (e.g., haemorrhages that exhibit the same imaging characteristics as vessels and confounds the automatic segmentation process) could be contributing to the lower Grad α observed. Indeed, Grad α was significantly lower in PDR when comparing to grade 0, 2 and 3  Table 3. The comparison of alpha-shape descriptors and FD obtained from manual annotation and automatic segmentation. DSC dice similarity coefficient with manual annotation as reference, MRE mean of the relative error with manual annotation as reference, Max maximum relative error with manual annotation as reference. n, number of images. + Pearson correlation test with null hypothesis that the correlation coefficient is zero. *Wilcoxon rank sum with null hypothesis that medians are equal between groups. www.nature.com/scientificreports/ but was not observed in grade 1. Individuals with grade 1 DR contain a small number of microaneurysms and is therefore expected to be vastly different from PDR. Grade 1 was the second smallest group in the MESSIDOR dataset and may be underpowered to detect a significant difference between PDR and grade 1. However, assuming that the dataset was sufficiently powered other reasons could include image quality such as the presence of artefact that may be contributing to a similar distribution of Gradα measured from PDR cases. To address such limitations, manual correction of the vessel segmentation or improvements to automatic vessel segmentation would be required. Nevertheless, Grad α could be a potential metric to evaluate segmentation quality in datasets that are several orders of magnitude greater than MESSIDOR where visual assessment and manual correction of the vessel skeleton prior to quantification of morphology is not feasible. We achieved good classification performance for PDR (AUC 0.67 to 0.76) that was improved with the addition of FD as a feature into the model and was comparable to the performance achieved by Orlando et al. 13 . Orlando et al. further improved their performance to AUC 0.91 when red lesion features were included as features for PDR classification 13 . We used fivefold cross validation and bootstrap to obtain performance measures and feature importance, when larger PDR datasets become available a held out experimental set-up should be implemented along with the addition of red lesion features if available.
We conducted a comparison of α-shape descriptors and FD measured from manual and automatic segmentations to assess measurement accuracy. The main differences between the manual and automatic segmentations are the presence or absence of small vessels. Additionally, manual annotations in the HRF are meticulously drawn, such that vessels are annotated to the complete length of the path. As vessels taper towards the end, they can lose contrast and may not be detected by the algorithm. In previous work that compared FD measured from automatic segmentations to manual annotations from two human observers (DRIVE dataset) reported mean relative errors between 7.1 and 9.3% and attributed this to the accuracy of the segmentation method 17 . We showed that α-shape parameters are also affected by segmentation accuracy. We therefore recommend careful selection of segmentation algorithms (e.g., use of multiple algorithms to assess reproducibility) to ensure α-shape and FD measures are as reliable as possible. Experimental results did not show a significant difference between healthy and DR that could be due to low sample size and presence of a broad range of DR severities in the DR group. FD measures obtained from manual and automatic segmentations showed a significant difference between groups, suggesting that it may be more sensitive than α-shape descriptors.
Interpretation of objective measures derived from the automatic segmentation of the retinal vasculature requires assessment of statistical significance of difference between groups and their clinical importance. To establish a threshold that identifies a change to be clinically important requires group-level between patient contrasts and within-patient longitudinal contrasts and therefore could not be established here. Nevertheless, providing careful segmentation procedures are followed, α-shapes appear to be a promising measure of vascular geometry. Future applications of α-shapes for describing retinal vasculature morphology could be applied to research involving the assessment of vascular changes in other ocular diseases as well as characterising vascular changes related to systemic comorbidities or neurodegenerative diseases that are increasingly reported to manifest in the retina 33,34 .

Data availability
The datasets generated during and/or analysed are available from the author on reasonable request.