Automated optic disk segmentation for optic disk edema classification using factorized gradient vector flow

One significant ocular symptom of neuro-ophthalmic disorders of the optic disk (OD) is optic disk edema (ODE). The etiologies of ODE are broad, with various symptoms and effects. Early detection of ODE can prevent potential vision loss and fatal vision problems. The texture of edematous OD significantly differs from the non-edematous OD in retinal images. As a result, techniques that usually work for non-edematous cases may not work well for edematous cases. We propose a fully automatic OD classification of edematous and non-edematous OD on fundus image collections containing a mixture of edematous and non-edematous ODs. The proposed algorithm involved localization, segmentation, and classification of edematous and non-edematous OD. The factorized gradient vector flow (FGVF) was used to segment the ODs. The OD type was classified using a linear support vector machine (SVM) based on 27 features extracted from the vessels, GLCM, color, and intensity line profile. The proposed method was tested on 295 images with 146 edematous cases and 149 non-edematous cases from three datasets. The segmentation achieves an average precision of 88.41%, recall of 89.35%, and F1-Score of 86.53%. The average classification accuracy is 99.40% and outperforms the state-of-the-art method by 3.43%.

www.nature.com/scientificreports/99.17%, and 99.08%.Another used approach is machine learning.Fatima et al. 6 developed a hybrid feature-based papilledema detection system.They first manually detected the OD region.The SVM classifier with thirteen extracted features from color, GLCM, statistical features, and intensity line profile was then used for classification.The method was evaluated on a small subset of the STARE dataset comprising 20 swelling cases and ten normal cases.The reported average values of the performance measures were 100% sensitivity, 95% specificity, 91.67% precision, and 96.67% accuracy, respectively.Yousaf et al. 7 extracted six related vascular features and four GLCM textual features from 36 manually cropped non-edematous OD boundaries and ten edematous OD boundaries from the STARE dataset.The classification was performed using the supervised support vector machine (SVM) classifier with a radial basic function kernel.They reported accuracy, sensitivity, precision, and specificity measurements of individual features of 95.65%, 100%, 83.30%, and 94.40%, respectively.In both studies in the machine learning approach, the OD region used as a domain for feature extraction was manually cropped.In addition, each experiment was done on a single small dataset with an imbalance of non-edematous and edematous ODs.
As automatic OD localization and segmentation are essential steps in our work, reviews of these tasks are also provided.Many localization techniques have been proposed based on the optic disk intensity, shape, size, color, and vessel information.We have summarized the work related to OD localization in Table 1.Although several techniques with many different features have been used for OD segmentation in the past, and some even achieved accuracy as high as one hundred percent, those methods are evaluated on collections in which most images are non-edematous OD.They did not consider edematous cases.All OD localization and segmentation methods rely on typical normal ODs' intensity, shape, and size.Therefore, these methods were not suitable for edematous ODs.Thus, the accuracy may not be as good as they claimed when involved with more ODE cases.This is because of physical changes in OD appearances, such as color, brightness, and size, and associating vessel structures that tend to be incomplete and tortuous.Reviews of automatic classification, localization, and segmentation of edematous OD are provided in the next section.The summary of the techniques used for OD segmentation is shown in Table 2.

Objectives, novelty, and contributions
This study extended our previous work 8 , initially presented at the 19th International Conference on Electrical Engineering, Computer, Telecommunications, and Information Technology (ECTI-CON 2022).In that conference work, we initially introduced a factorized gradient vector flow (FGVF) 9 , a special kind of gradient vector flow for texture segmentation, to segment the edematous ODs.It is experimentally proven on a small dataset containing 35 images to yield high performance.
The extension parts from the previous work can be summarized as follows.
1. We experimented with using FGVF to segment the OD for non-edematous ODs.The previous work was only done on edematous ODs. 2. We experimented with more images from two additional public datasets containing both types of ODs.The total number of images used in the experiment is 295 with 146 edematous and 149 non-edematous cases.3. We demonstrated that the precise OD boundary is useful for edematous OD classification.
The use of FGVF to identify the boundary of the optic disk (OD) is a groundbreaking technique in the field of research.This approach shows great promise for advancing OD detection methodologies, especially in cases where ODs are swollen.In comparison to four other state-of-the-art methods, experiments show that FGVF can provide precise OD segmentation results, regardless of whether the OD is edematous or non-edematous.This is a significant finding in the field of ophthalmic image processing, where images with mixed types of ODs are common.Moreover, the classification of OD types can be particularly useful for ODE prescreening.

OD localization
To locate the OD in cases where the vascular networks were incomplete, we used the hybrid localization method (HLM) 24 .This was because the vascular networks were often incomplete in edematous cases.The HLM method was effective on all types of vascular networks, regardless of their structural completeness.The HLM method first analyzed the structure of the vascular network.If the vascular network was complete, the main vessels appeared in a horizontal parabolic shape.The HLM method assumed that the vertex of the parabola was the location of the OD.If the vascular network was incomplete, it appeared as several broken lines.The OD location is determined

OD boundary segmentation
The OD region of interest (ROI) was first defined as a square centered at the OD location obtained from the HLM method.According to the size of OD in our datasets, we assumed that the diameter of the non-edematous OD was one-sixth of the retina's diameter 24 .As the edematous OD's region was commonly larger than the average size of the non-edematous OD, the ROI square's width of both edematous and non-edematous OD was set to one-third of the retina's diameter.Figure 4a and b depicts the original image and the ROI region.Next, the image contrast was enhanced.The color space transform was applied to convert from an RGB to a L*a*b* channel.The Contrast-Limited Adaptive Histogram Equalization (CLAHE) 8 was performed on the L channel on the ROI.For removing vessels from ROI, a masked image was created by multiplying the green channel of the original image with the binary ROI.Gaussian filtering was applied to the green channel of the masked image to smooth it.Then, the region filling was performed on the pixels within the mask based on Laplace's equation, removing vessel-related regions.Figure 4c and d illustrate the results after these processes.
A factorized Gradient Vector Flow (FGVF) proposed by Gao et al. 9 was employed in our work to segment the OD boundary.The following FGVF pseudocode illustrates the main tasks in a recursive manner.The initial contour (C) was defined as a circle with a radius of 1/4 of the retina's width, centered at the OD location.Subsequently, the texture feature matrix Y was computed from the vessel-removed image from the prior step.The FGVF algorithm takes the texture feature matrix Y, the initial contour C, the number of rounds i as inputs.It repeatedly evolves the contour C to be closer to the OD boundary until convergence.
The FGVF algorithm uses the following functions.MakeFeatureMatrix(Img) takes the image Img as an input.It returned a texture feature matrix calculated by using local spectral histograms 39 and a factorized-based texture segmentation method proposed by Yuan et al. 40 .
PerformEvolution(Y, C) takes a texture feature matrix Y and a contour C. It evolves using the level set function and returns the new contour.The contour evolution is performed using level-set regularization proposed by Li et al. 41 .
CheckConvergence(C, C*) takes contours C and C* as inputs.If the average differences along the x and y directions between the input contours are less than a convergence threshold, the function returns true; otherwise, it returns false.In our experiment, we used 0.05 for a convergence threshold.
Figure 4e-h displays contours at the initial round, 100th round, 400th round, and upon convergence.
The contour evolution based on the factorization-based fitting energy and level-set regularization can be mathematically explained.Given φ a signed distance function of a contour curve and R is the presentation feature 9 .
The FGVF energy function (E FGVF ) consists of two energy terms: a factorization-based fitting energy ( E data ) 9 and a level-set regularization term ( E regularization ) 41 .
where τ and υ are two positive constants to control the proportion of E data and E regularization .In our edematous and non-edematous OD boundary segmentation, we set the constant values τ = 50 and υ = 1.5.These values are tested empirically to yield the best result.
The first energy term E data is derived from the matrix factorization techniques.Equations ( 2)-( 5) collectively contribute to E data .Terms A and B are determined using the Heaviside function and the weight vectors ω o and ω b .The weights are calculated from the presentation feature R and the feature matrix Y.
(1) www.nature.com/scientificreports/where Ω is a 2D image domain, x is a point in the domain, Ω o and Ω b are defined as the object region and the background region (i.e.Ω = Ω o ∪ Ω b ) , H ǫ (φ) is a Heaviside function, ω o and ω b are the weights of the object and background regions, R is the representative features, Y is the feature matrix of the ROI region calculated using factorization based method for textual image segmentation proposed by Yaun et al. 40 , β is a matrix whose columns are region weight vectors, and ǫ is the additive noise.In our work, ǫ is set to 0.5.The object and background of ROI are divided into two parts with different textural feature maps.The Y in Eq. ( 6) refers to the resultant matrix of the MakeFeatureMatrix(Img) function in the prior FGVF pseudocode.
The second energy term E regularization shown in Eq. ( 1) is expressed as: where ∇ φ is the derivate of the level set function, the deformation process is repeated until the contour converges into the object boundary.Equation ( 7) enforces regularization and smoothness in the level-set function.
The update of the level set function at each round is described in Eq. ( 8).
The evolving contour is a level set of φ , expressed as in Eq. ( 9).
This C in Eq. 9 refers to the resultant contour of the PerformEvolution (Y, C) in FGVF pseudocode.

Edematous classification
A compilation of 40 different appearance-based and statistical-based features of OD was made from various literature sources.The maximum relevance minimum redundancy (mRMR) algorithm 42 was then utilized to pick the most relevant and non-redundant features from the initial list.The mRMR algorithm prioritized features that offer informative data while minimizing any redundant information.It finally picked a set of 27 features, which we grouped into four categories, namely GLCM, vessel, color, and intensity line profiles.Below is a detailed list of the selected features.

Gray-level co-occurrence matrix features
Gray-level co-occurrence matrix (GLCM) is a statistical technique for analyzing texture that considers the spatial relationship of pixels 43,44 .GLCM calculates the texture based on pairs of pixels with specific values and their spatial arrangement.Ten GLCM features are extracted.Let M be a co-occurrence matrix with N dimension, where N is the number of gray-values, all pairs of intensities i, j are its coefficients and coordinates of the elements, p is the normalized co-occurrence matrix, µ x , µ y and σ x , σ y are the mean and standard deviations for the matrix p's rows and columns, respectively.
1. Autocorrelation (autoc) computed as the sum of the product of each element in the matrix p and the product of their distance from the mean that refer to the absolute differences between the row and column indices of the element in p and the mean row and column indices, respectively.It is high in edematous OD due to having similar intensity values, while non-edematous OD has lower autocorrelation because of high intensity change between optic disc and optic cup of the normal condition.
2. Contrast (contr) measures the difference in color shades and brightness of the region.A higher contrast value indicates a higher variation in gray level between neighboring pixels.Thus, non-edematous OD has larger contrast value than the edematous condition.
where n = i − j and N g is quantized gray levels.

Correlation (corrp) uses means and standard deviations to quantify the linear relationship between pixel
intensities in the matrix p .The low variations in pixel intensities of the edematous case show high correla- tion.
4. Cluster prominence (cprom) measures the presence of clusters in the image, where a higher value indicates a greater prominence of clusters in the image.Thus, non-edematous OD has high value and edematous OD has low value.www.nature.com/scientificreports/ 5. Cluster shade (cshad) measures the degree of asymmetry in the grayscale pair distribution.Non-edematous condition has high asymmetry in the distribution of the matrix p and edematous case perform low asym- metry.
6. Dissimilarity (dissi) measures the average absolute differences between pixel intensities in the matrix p .When non-edematous OD has significant changes in texture of optic cup and disc, the value is high.
7. Energy (energy) measures the uniformity of the image pixels.Texture is likely uniform in edematous OD and varying in non-edematous OD.Thus, energy value increase in abnormal disruptions of texture patterns.
8. Entropy (entro) measures the disorder in the distribution of pixel pairs.The value is high when the matrix p 's elements are uniformly distributed.The homogenous characteristic of edematous OD condition has low entropy.9. Homogeneity (homop) measures image homogeneity with larger values for smaller gray tone differences in pair object.The non-edematous OD has low homogeneity compared to the edematous OD.
10. Max probability (maxpr) measures the most frequently occurring intensity pair in the image.The edematous ODs tend to have a lower value of maxpr than the nonedematous.

Vessel features
The following are the vessel features used in the experiment.
1. Vessel disk continuity Index (VDI) is the number of disjointed vessel regions in the segmented vascular network of OD images.Non-edematous OD image usually has a completely connected vascular structure, resulting in a low VDI value.In comparison, an edematous OD image usually has more broken vessels, especially in a severe case, resulting in a higher VDI value 3 .2. Vessel disk continuity index to disk proximity (VDIP) is a VDI that calculates within the scope of the segmented OD region.3.An area of the largest vessel region is the number of pixels in the largest vessel region.It offers details regarding the prevalence or range of the largest vessel structure.The edematous ODs tend to have a smaller number than the non-edematous ODs due to less completeness of the vascular network.4. Mean vessel area-The ratio of the sum of vessel pixels to the total number of connected vessel clusters.
Edematous ODs usually have this number lower than non-edematous ODs due to the vessel compression effect.5.A standard deviation ( σ ) of the probability of intensity distribution is defined as follows.
where N is the total number of pixels in the image, x represents each pixel intensity value, µ is the mean of image intensity distribution.The σ values of the edematous ODs tend to be higher than the non-edematous ODs. 6. Kurtosis distribution ( κ ) is a measure of the tailedness of an intensity distribution defined as follows.Generally, the averages, the maximums, and the standard deviations of the intensity of the non-edematous ODs are larger than the edematous ODs.In contrast, the minimum intensity values of non-edematous ODs are lower than those of edematous ODs.

Datasets, classifier, and evaluation
The programs were implemented using MATLAB R2022a and ran on DELL IN5406 (Intel Core i7-1165G7 Processor).The experiments were tested on three datasets.The first dataset downloaded images from the Internet 45,46 includes 35 edematous and 38 non-edematous ODs images with the dimensions between 600 × 600 and 2300 × 1900.The selected fundus images with optic disk edema from RFMiD public dataset 47 contained 91 edematous and 91 non-edematous ODs images with dimensions between 2144 × 1424 and 4288 × 2848.From the RFMiD2.0public dataset 48 , 20 edematous and 20 non-edematous OD images with the dimensions 512 × 512 and 2048 × 1536 were selected.A total of 295 OD images with 146 edematous and 149 non-edematous cases were used in the experiments.
For the ODE classification, we selected a Linear Support Vector Machine (SVM) since it is effective with datasets with many features.To minimize over-fitting, we used fivefold cross-validation approach with 80% training and 20% testing.However, when dealing with a new image with different characteristics from the current dataset, over-fitting may still occur.Additionally, it is important to note that there are limited publicly available retinal images with edematous ODs.Thus, it is currently not possible to solve the issue of over-fitting by simply increasing the size of the dataset.
For OD localization, the performance was measured using a location accuracy (Acc loc ) defined in Eq. ( 22).
where C is the number of images the method correctly localizes the OD, and N is the number of images.Remark that the case is successful when the method's calculated OD location is within the ground truth contour.The performance of the OD segmentation method was evaluated using precision, recall, and F1 measures.The evaluation formulas are shown in Eqs. ( 23)- (25).
where TP, FP, TN, and FN are the number of pixels that are true positive, false positive, true negative, and false negative, respectively.
For edematous classification, we compared the performances of each feature and all together features using a support vector machine (SVM) linear classifier.The classification accuracy (Acc classify ) is defined in Eq. ( 26).
where C Ede and C Non are the numbers of images correctly classified as edematous and non-edematous and N is the number of images.

Numerical results and discussion
This section presents comparative and quantitative studies of localization, segmentation, and classification of edematous OD compared to the existing methods.

OD localization
We compared the hybrid localization method (HLM) 24 used by our method against the feature projection (FP) 22 and adaptive thresholding (AT) 10 methods.Selected cases of localization results from non-edematous and edematous groups of two datasets are depicted in Fig. 6.
The numerical results are reported in Table 3.For non-edematous ODs, most methods could locate the OD efficiently.The FP performed worse than other methods because it relied only on vessels.When the vascular network was incomplete in some edematous cases, the FP failed.AT could sometimes spot abnormally high bright spots as the OD.
Results from Fig. 6 and Table 3 showed that the HLM method used by our algorithm achieved the best average Acc loc of 97.88% for all three datasets and was considerably higher than the FP and AT methods.The Acc loc values of all three methods were lower in the edematous cases than in the non-edematous cases.Across all OD types, the average Acc loc of HLM was higher than FP and AT by 12.48% and 6.04%, respectively.Moreover, HLM localization performance was significantly superior to the other two comparative methods, especially in edematous cases.For such cases, Acc loc of HLM was higher than FP and AT by as much as 22.49% and 12.03%, respectively.( 22)

OD segmentation
We compared the factorized gradient vector flow (FGVF) 8,9 used in our work against four other comparative methods: alternated deflation-inflation gradient vector flow (ADI-GVF) 37 , traditional gradient vector flow (GVF) 36 , region growing (RG) 29 , and super-pixel clustering (SPC) 30 .All methods except super-pixel clustering required initial points.The OD locations obtained from the HLM method were the initial points.Figure 7 shows examples of segmentation results from different approaches for edematous and non-edematous ODs.
Most methods performed better on the non-edematous ODs than the edematous ODs.For edematous OD, the methods in the GVF family showed undersegmentation, while the region growing and superpixel clustering showed oversegmentation.Most methods worked well for non-edematous OD.Table 4 shows the numerical performance comparison of segmentation methods.
The following findings can be summarized from the results of Table 4.
In the case of non-edematous images, both GVF and FGVF methods have F1 measures that are significantly higher than other comparative methods.On average, the improvement of FGVF over the second-best method (GVF) is only 0.21%.However, FGVF outperforms the poorest method (ADI-GVF) by 21.16%.
It was found that for images with edema, all methods performed worse than those without edema.Among all the methods, FGVF was the best and had significantly better results than GVF and other methods.On average,  www.nature.com/scientificreports/ the improvement of FGVF over the second-best method (GVF) was 2.99%, while the improvement of FGVF over the poorest method (RG) was 15.54%.
In general, regarding mix cases, both GVF and FGVF have F1 measure values that are fairly close, but significantly better than other methods.Precision-wise, GVF was slightly better than FGVF, but FGVF had considerably better recall than GVF.This resulted in FGVF having a better overall F1 measure than GVF.However, the ADI-GVF method was the poorest performer among them.
The RFMiD2.0 dataset is known to be more challenging for most methods due to the low resolution and indistinct OD region in edematous OD images.However, when only considering edematous OD images in this dataset, FGVF still performs best.

Edema classification
Table 5 compares the linear SVM classifier accuracy of classification performance using each sole feature set and combined feature sets on different datasets.
The findings and discussions from Table 5 are as follows.
1.The proposed method achieved an average accuracy of 99.40%, which was the highest accuracy recorded.
2. The results of the average accuracy classification for a single feature type showed that the intensity line profile, vessel, color, and GLCM gave the best to worst results, respectively.However, in the proposed work, the average accuracy significantly improved.The overall improvement of all feature types combined compared to the best type (the intensity line profile) was 2.6%, while compared to the worst type (GLCM) was 11.4%.3. It should be noted that our proposed method had shown better performance than the method proposed by Yousaf et al. 7 by 3.34%.This improvement could be attributed to the fact that Yousaf et al. used only ten features from vessels and GLCM, while our method also employs features from the Color and intensity line profiles.This suggests that using additional features could help improve the classification results.4.After analyzing the unsuccessful cases, we found that the accuracy of classification depended on several factors, such as the stages of edema in the dataset and the appearance of the OD.We noticed that when the images dealt with the mild edema stage in the dataset, the classification accuracy was lower.This was because the differences in characteristics between mild edema and normal OD showed minimal changes in the appearance of the disk.Additionally, some non-edematous OD with unclear boundary resulted in incorrect segmentation of the OD, which led to extracting wrong features and consequently resulted in inaccurate classification.Figure 8 shows examples of an edematous OD misclassified as non-edematous (false negative) and a non-edematous image misclassified as edematous (false positive).

Conclusion
This paper presents an automatic classification and segmentation of optic disks with edematous and non-edematous based on the FGVF segmentation model using HLM initialization and classification results from a linear SVM classifier.The proposed method was evaluated on 146 edematous and 149 non-edematous images from Internet and RFMiD datasets by comparing the proposed localization, segmentation, and classification performances against the existing methods.The HLM worked well for OD localization and correctly located the OD in 295 out of 292 images with 97.88% accuracy.The proposed FGVF achieved an average segmentation precision of 86.56%, recall of 88.19%, and F1-score of 86.48%.The average classification accuracy was 99.40%.However, the FGVF method used in the OD segmentation algorithm had limitations, including high computational demands  www.nature.com/scientificreports/and sensitivity to initial conditions.For edematous OD classification, accuracy relied heavily on the precision of OD segmentation.Finding more useful features, such as the cloud OD boundary and the ratio of OD diameter to that of the retina, and improving the limitations of FGVF will be our future work.

Figure 3 .
Figure 3. Illustration of the HLM method used for OD localization (rectangle) in non-edematous (left) and edematous ODs (right).

Figure 4 .
Figure 4. FGVF procedure illustration (a) Original Image, (b) Region of Interest (c) Contrast-Enhanced Image, (d) Vessel Removed Image, (e) Seed point and initial contour (f) after the 100th round of FGVF's evolution process, (g) after the 400th round of FGVF's evolution process, (h) OD boundary after FGVF's convergence.

Figure 6 .
Figure 6.Examples of OD localization results of edematous cases from FP (yellow square), AT (red circle), and HLM used by the FGVF method (blue hexagram) for non-edematous (top) and edematous (bottom).

Figure 8 .
Figure 8. Examples of false negative (left) and false positive (right) cases.The black solid contour is the ground truth, the hexagon is the OD location, and the dash line is the OD boundary.

Table 1 .
Reviews on OD localization techniques.

Table 2 .
Reviews of OD segmentation techniques.Fang et al. 26 Using biregional contour evolution model from the two-level set functions.The intensity, edge, and area features are considered in the method, and the Edge indicator function (EIF) is computed to differentiate OD and OC edges 1341 images from Dhristi-GS, DRIVE, REFUGE Jaccard 93.20, Dice 96.48, Accuracy 99.77 34Using a combination of the three energies: phasebased boundary, PCA-based shape, and region energies 1409 images from MESSIDOR, ONHSD, DRIONS Overlap 90.54 Xue et al.33The hybrid level set model (HLSM) included distance-regularized, line integral and area integral, area-based, and shape-based models 138 images from DRSHTI-GS, TMUEH Intersection over union 92.75, Four-side evaluation 464.36   Gao et al.34Using saliency detection and thresholding techniques to get a rough OD boundary. Thval fitting model is used to segment higher accurate boundary 229 images from DIARETDB0, DRSHTI-GS Overlap 66.59, Accuracy 96.30, F1-score 95.1 Abdullah et al. 35 Using the fuzzy clustering mean method to localize the location.The active contour model is Vol.:(0123456789)Scientific Reports | (2024) 14:371 | https://doi.org/10.1038/s41598-023-50908-5

Table 3 .
Comparison of the OD localization performance of FP, AT, and HLM.The highest number in class is bold.

Table 4 .
OD segmentation performance.*The highest is bold.

Table 5 .
7ccuracy comparisons of the proposed method (all featured combined) against each feature set and also against a state-of-the-art method (Yousaf et al.7).*The highest number in class is bold.