Robust Detection and Modeling of the Major Temporal Arcade in Retinal Fundus Images

: The Major Temporal Arcade (MTA) is a critical component of the retinal structure that facilitates clinical diagnosis and monitoring of various ocular pathologies. Although recent works have addressed the quantitative analysis of the MTA through parametric modeling, their efforts are strongly based on an assumption of symmetry in the MTA shape. This work presents a robust method for the detection and piecewise parametric modeling of the MTA in fundus images. The model consists of a piecewise parametric curve with the ability to consider both symmetric and asymmetric scenarios. In an initial stage, multiple models are built from random blood vessel points taken from the blood-vessel segmented retinal image, following a weighted-RANSAC strategy. To choose the ﬁnal model, the algorithm extracts blood-vessel width and grayscale-intensity features and merges them to obtain a coarse MTA probability function, which is used to weight the percentage of inlier points for each model. This procedure promotes selecting a model based on points with high MTA probability. Experimental results in the public benchmark dataset Digital Retinal Images for Vessel Extraction (DRIVE), for which manual MTA delineations have been prepared, indicate that the proposed method outperforms existing approaches with a balanced Accuracy of 0.7067, Mean Distance to Closest Point of 7.40 pixels, and Hausdorff Distance of 27.96 pixels, while demonstrating competitive results in terms of execution time (9.93 s per image).


Introduction
The analysis of vascular structures in the retina can facilitate the monitoring and diagnosis of different types of ocular pathologies. The Major Temporal Arcade (MTA) is the thickest vascular structure in the retina and it is composed of the superior and inferior temporal arcades [1]. In clinical practice, the visual examination of the morphological integrity of the MTA and the angle between the two temporal arcades, also called Temporal Arcade Angle (TAA), are used as an indicator of the severity of diabetic retinopathy, myopia, or hypertension. Consequently, the numerical modeling of the MTA plays an essential role in achieving its quantitative analysis, while improving medical diagnosis.
In literature, a handful of works have addressed the analysis of the MTA. Oloumi et al. [2] introduced image processing techniques for the segmentation, tracking, and measurement of the MTA width in order to detect plus disease and retinopathy of prematurity. Nabi et al. [3] proposed a two-step method for the identification of the MTA. Firstly, the retinal vessels are detected through Gabor filters. Then, the Hough transform along with graph theory were applied to separate the MTA. Fleming et al. [4,5] considered the MTA length as a criterion to measure the quality of fundus images, along with other elements such as the optic disc location and the macula size. The MTA was automatically determined by semi-elliptical templates in a range of sizes, using the generalized Hough transform. Fledelius and Goldschmidt [6] studied the changes in the MTA geometry for patients with high myopia, finding a correlation between the increase in myopia and the reduction in the TAA. The analysis was performed manually by ophthalmologists, who located MTA and measured the TAA.
To simplify the task of automatic tracking of MTA morphological alterations over time, Oloumi et al. [7] proposed the parameterization of the MTA by using a parabolic model. A weighted version of the Hough Transform was applied to fit a parabola on the detected vascular tree. A similar approach was presented by Oloumi et al. [8], where a dual parabolic model was designed (i.e., two parabolas adjusted to the upper and lower parts of the MTA) as a solution to the problem of MTA non-symmetry for some fundus images, especially those belonging to sick patients. These models were used to prove that the openness of the MTA decreases when patients suffer from Proliferative Diabetic Retinopathy and Retinopathy of Prematurity [9].
More recently, population-based methods have been presented as an alternative for the detection of parabolic objects in medical images. Guerrero-Turrubiates et al. [10] introduced a method that applied a Univariate Marginal Distribution Algorithm (UMDA) to approximate a parabola for retinal images. The approach creates individuals by concatenating three pixel indices in the image domain. Then, each individual is evaluated using as a fitness function the Hadamard product between the input image and the resulting parabola image. Valdez et al. [11] developed a parabola detection algorithm for the MTA localization in fundus images. Instead of using the Hough transform, which is computationally expensive, a fast hybrid method was proposed that combined the UMDA algorithm with the Simulated Annealing (SA) strategy, which guided the search towards promising regions. The objective function used a segmented image of the vascular structure, weighing pixels according to the distance to the parabola vertex. These two UMDA-based methods obtained superior performance in terms of computational time compared to Hough-based methods.
In this paper, a robust method for the detection and piecewise parametric modeling of the MTA in fundus images is presented. The algorithm follows a weighted-RANSAC strategy, building multiple MTA models from random blood-vessel points taken from the blood-vessel segmented retinal image. Each model consists of a piecewise parametric curve with the ability to consider both symmetric and asymmetric scenarios. To choose the best MTA model, the method considers blood-vessel-width and foreground-location features, extracted through the Distance Transform and grayscale intensities, respectively. Both attributes are merged and used as a probability function to weight the percentage of inlier points for each model. This procedure promotes selecting a model based on points with high MTA probability.
The contributions of this work are summarized as follows: 1.
Blood-vessel-width and intensity features are integrated into the MTA detection and modeling process.

2.
A modeling strategy addressing both symmetric and asymmetric scenarios is presented to improve the MTA characterization.

3.
A weighted-RANSAC scheme is included to robustly select a MTA model configuration.

4.
A set of MTA manual delineations for the benchmark DRIVE dataset has been released for scientific purposes.
The rest of this paper is organized as follows. In Section 2, the proposed method is explained in detail. The dataset, evaluation metrics and experimental results are discussed in Section 3. Finally, in Section 4, the conclusions obtained from this work are presented.

Methods
The proposed method consists of the following steps: (1) pre-processing of the raw images to deal with the illumination and contrast changes expected in medical images; (2) automatic segmentation of the vascular tree to reduce the MTA search space; (3) the extraction of features from the segmentation carried out in the previous step; and (4) the construction of a numerical model of the MTA from the extracted features. Each of these steps is described in detail in this section.

Pre-Processing
Fundus images are affected by uneven illumination and low contrast, producing variations in pixel intensity [12][13][14]. In this work, the original RGB fundus image is converted to grayscale and then pre-processed through the Contrast Limited Adaptive Histogram Equalization (CLAHE) algorithm, which addresses these intensity changes and improves blood-vessel segmentation accuracy by applying a re-distribution of their intensity histogram [15][16][17].
The CLAHE algorithm partitions the image into small rectangular regions and applies a histogram-based local-contrast enhancement procedure. First, the intensity occurrences that exceed a particular value in the histogram of each region in the image are truncated. Then, the truncated occurrences are redistributed uniformly over the histogram. This redistribution of occurrences results in an increase in contrast between the background and the objects of interest, while the truncation procedure aims to avoid noise amplification problems, which are typical issues of traditional histogram equalization methods. Further information on this pre-processing technique for retinal images can be found in [18][19][20].

Automatic Segmentation of the Vascular Tree
The inclusion of the automatic blood-vessel segmentation step in the proposed algorithm aims to delimit the MTA search space efficiently. Hence, a good trade-off between high accuracy and computational time must be considered in order to select the adequate approach for this step.
In literature, the problem has been widely addressed by unsupervised methods [21][22][23][24], machine learning strategies [25][26][27], and deep learning techniques [28][29][30][31]. Many recent models improve blood-vessel segmentation. However, their advances are generally correlated to a greater number of parameters or increased numerical complexity. Moreover, these methods focus on refining specific cases, such as thin vessel detection, which is less relevant in the MTA detection task.
Considering this context, a Convolutional Neural Network (CNN) with the U-Net architecture has been adopted to perform the task. This approach has proven high efficiency for biomedical image segmentation on reasonable computational time [32][33][34][35]. Moreover, this model constitutes an end-to-end system with learnable parameters and few hyper-parameters to calibrate. Unlike the original four-level model with an initial depth of 64 channels, a simplified three-level model with an initial depth of 32 channels is employed here. The architecture is assembled by the contracting path (or encoder), the latent path (or bottleneck), and the expanding path (or decoder). Skip-connections are added between feature maps from encoder to decoder to propagate previous information to deeper network layers. Figure 1 illustrates the overall design. Each level in the contracting path consists of a double 3 × 3 convolutional layer + ReLu block and a Maxpooling layer. Similarly, each level in the expanding path consists of a double 3 × 3 convolutional layer + ReLu block followed by a Upsampling layer. A dropout layer has been added between each block to improve generalization [36]. Moreover, a patch-based approach has been adopted in the design to increase the number of images available for training: random 48 × 48 pixel patches are extracted from the grayscale pre-processed train images and are used as input for training the model. Figure 2 shows an overview of this process. For the testing part, ordered 48 × 48 pixel patches with an overlapping of five pixels are extracted from the test images. The final prediction is obtained by averaging the predictions made over each pixel.
Binary cross entropy has been selected as a loss function to adjust the network parameters.
where θ represents the parameters of the network, m is the number of samples or pixels to classify, y i is the true label of sample i andŷ i is its predicted label, i.e., the network output.

Feature Extraction
The MTA is the thickest blood vessel that appears in the foreground of the fundus image. The presented method proposes to quantify blood-vessel thickness and its location in the foreground of the image, so that both attributes contribute to a robust MTA detection and modeling.

Vessel Thickness
Determining blood-vessel-width is a difficult task in fundus retinal imaging due to the variety of blood-vessel amplitudes. A quick and indirect measurement can be performed by using the Distance Transform [37][38][39], since it results in a map containing each blood-vessel pixel's distance to its nearest background pixel.
The two-dimensional distance transform can be formulated as follows. Let L be the set of sites or pixels in a two-dimensional binary image I ∈ R N×M , and let A and B be two non-overlapping subsets of L, i.e., L = A ∪ B and A ∩ B = ∅, such that set A contains all the sites u = (x, y), u ∈ L that belong to the foreground of the image and set B contains all the sites u ∈ L that belong to the background of the image: The Distance Transform is a function that generates a map D ∈ R N×M for which each value at site u corresponds to the smallest distance from u to B: where d(·, ·) is a distance metric between two points, usually the Euclidean distance metric. Through a normalization process, i.e., dividing each element in D by the maximum value, and using the segmented image of the vessels as input, the Distance Transform can be used as a metric to determine the vessel width in each position of the image, allowing for the distinction of the following cases: • When a pixel is not located in a blood vessel, the metric returns a value of zero. • When a pixel is located in the center of a thin blood vessel, the metric returns a value close to zero. • When a pixel is located in the center of a thick blood vessel, the metric returns a value close to one.
The metric has a drawback in determining blood-vessel width: as the pixels get closer to the edges, their value gets closer to zero, even when they belong to thick blood vessels.
A strategy to avoid these misleading values is to preserve only the center-line pixels of the vascular structure. However, this approach alters the percentage of pixels belonging to thick blood vessels in the image, which may lead to performance decreasing in the numerical modeling procedure.

Foreground Location
Notice that having an indicator of blood-vessel width is helpful to some extent in locating sites of the image that may belong to the MTA. However, considering that some thick blood vessels appearing in the background of the image could also obtain high responses in the Distance Transform, additional information is required to make a more robust MTA detection.
The blood vessels in the foreground can be identified because their intensities are visibly darker than those that appear in a second plane. Following this idea, a naive foreground location map F ∈ R N×M can be calculated by taking the complement of the normalized grayscale blood-vessel intensity image G ∈ R N×M .
In the resulting map, the highest values correspond to pixels belonging to firstplane blood vessels, while lower values correspond to pixels belonging to background blood vessels.

MTA Probability Map
The foreground location and blood-vessel thickness feature maps are averaged in order to consider a joint contribution. The average image can be interpreted as a coarse MTA Probability Map P ∈ R, where each pixel represents the probability of that position in the image belonging to the structure of the MTA. The process of obtaining P is illustrated in Figure 3. Previous works regarding the MTA curvature approximation have presented models using parabolic Hough-based techniques that assume a symmetric shape. However, the MTA is not strictly symmetrical. It may present shape alterations derived from a disease or vascular damage. To address this issue, a piecewise parametric approach using quadratic spline curves is proposed.
The spline curve is a piecewise low-order polynomial function that approximates intrinsic shapes while avoiding abrupt oscillations (Runge's phenomenon) [40][41][42]. Consider a set of n + 1 ordered points or knots x 0 , x 1 , ..., x n and an integer k > 0 that have been specified. Let S(x) be a piecewise polynomial function defined on the interval [x 0 , x n ] as follows: where each piece p i with 0 ≤ i ≤ n − 1 is a polynomial function of degree at most k. To guarantee the continuity and smoothness of S(x), the two polynomials p i−1 and p i must share the values of their derivatives from the order (i.e., the value of the function) up to the derivative of order m at knot x i . Then, S(x) is said to be a spline curve of degree k and smoothness C m , or S ∈ C m in the neighborhood of x i . For instance, to build a model with a quadratic spline curve S ∈ C 1 from a set of ordered pairs u 0 , u 1 , ..., u n with u i = (x i , y i ) and 0 ≤ i ≤ n, a function S(x) as described in (5) must be defined, where the polynomial functions p i are of the form: Then, S(x) interpolates the given set of points, that is: The polynomials p i and p i+1 must interpolate the same value at the point x i+1 to ensure continuity along the interval [x 0 , x n ], thereby the following expression have to be satisfied: Furthermore, S (x) is also required to be continuous, i.e., The expressions in (7)- (9) can be used to construct an equation system to find the values of the 3n unknown coefficients of the polynomial functions p i , adding an additional constraint to the value of p 0 , typically p 0 = 0.

Weighted RANSAC
In the proposed method, the knots for the quadratic spline curve are randomly taken from the set of non-zero pixels in the MTA Probability Map obtained from the previous step. However, the MTA Probability Map values are noisy, since they were obtained from an estimation based on features extracted from the segmented image. A weighted RANSAC methodology is proposed to robustly select the points that best accomplish the particularities of the MTA.
The RANSAC method is a non-deterministic iterative algorithm for computing the parameters of a model M given a set of N points or observations that contain noise [43,44]. It consists of taking a minimum subset of points Q = u 1 , u 2 , ..., u h , h < N, necessary to compute the model parameters and observe the number of points that are well explained by it, i.e., the number of points that are within a margin distance from the points estimated by the model (inliers). This process is repeated a number of iterations or until a percentage of inliers is reached, keeping the model that best explains the dataset.
In the original RANSAC, each point inside the margin distance makes the same contribution to the inlier counting. In contrast, in a weighted-RANSAC scheme, the contribution w(u) of each point u = (x, y), u ∈ L to the counting is weighted using a criteria [45]. In the proposed method, the criteria is given by the MTA Probability Map P: where d(·, ·) represents a distance metric between two points -in this case, the distance between point u and the predictionû obtained with model Mand > 0 is the tolerance error value that determines if point u is considered as an inlier. Algorithm 1 describes in detail the weighted-RANSAC scheme. This strategy prefers models containing points at the center of thick vessels with strong intensities are preferred. This behavior is also imposed from the task of selecting points to estimate the model: although a uniform probability distribution is used to choose the points acting as knots, the proportion of points belonging to thick vessels is greater than the one corresponding to thin points. Hence, points belonging to thick vessels are selected more frequently.

Algorithm 1: The weighted-RANSAC algorithm
Data: Blood-vessel segmented image (I); MTA Probability Map (P ); Number of points(n); Polynomial degree (k); Margin distance ( ); Maximum number of iterations (max_iter); Result: Piecewise parametric model M; Model points model_points 1 max_inliers ← 0; 2 model_points ← ∅; 3 for i ← 1 to max_iter do 4 Create a set C of randomly select n blood-vessel points from I: Sort elements in C with respect to y-axis; 6 Compute model S with polynomial fuctions of degree k with respect to y-axis, using Equations (7)-(9) and set C; 7 Compute inlier count W as: with w(u) defined as in Equation (10)

Constraints on Point Selection
Point selection can result in models that do not cover the entire length of the MTA. On that account, a restriction has been added to choose points from three image regions (top, middle and bottom) in balanced parts, as shown in Figure 4. Under these conditions, the search for the optimal model with RANSAC is always conducted considering points around a larger area in the image, which improves the results of the MTA modeling. . Each point has the inlier weight that corresponds to its value on map P, represented here by its grayscale intensity.

Results and Discussion
The proposed method was evaluated through multiple comparisons with state-ofthe-art algorithms using the dataset DRIVE. First, the dataset DRIVE and its MTA manual delineations are presented. Secondly, the evaluation metrics are described. Third, the implementation details are explained. Finally, the comparative analysis is discussed.

Dataset and Delineation of the MTA
The dataset DRIVE [46], consisting of 40 retinal fundus images with size 565 × 584 pixels, is used to evaluate the performance of the proposed method. The partitions used for training and testing were taken as recommended by the dataset authors.
The DRIVE dataset has been designed for the blood-vessel segmentation task and does not provide ground-truth images for the MTA detection problem. Hence, hand-labeled MTA annotations have been created for this work by an expert ophthalmologist of the Ophthalmology Department of the Mexican Social Security Institute (IMSS) T1-León [47]. The set of MTA manual delineations is freely available (http://personal.cimat.mx:8181 /~ivan.cruz/Journals/MTA_drive.html, accessed on 11 April 2022). To the authors' best knowledge, this dataset is the first in the literature that releases MTA manual delineations for scientific purposes.

Evaluation Metrics
The metrics considered to evaluate the closeness between the model and the MTA ground-truth are Mean Distance to Closest Point (MDCP) and Hausdorff Distance, as proposed by Oloumi et al. [7].
MDCP measures the average of the distances of each point belonging to one set to the nearest point of the other set. Let A and B be two sets of points, where A is the set of points estimated by the model and B is the set of points in the ground-truth delineation. The MDCP can be defined as follows: where N is the cardinality of set A, a i its i-th element and DCP is the distance to the closest point, which is computed as follows: for j = 1, 2, ..., M, M being the cardinality of set B, b j its j-th element and || · || any norm operator, typically the Euclidean norm. Moreover, Hausdorff Distance also uses the previous definition of DCP to find the smallest distance of each point of A to the ground-truth set B, however it does not calculate an average. Instead, Hausdorff Distance takes the maximum DCP distance: In both measures, small values indicate that the model is a good fit for the ground-truth. The metrics regarding the MTA detection, taking into consideration that the model corresponds to a one-pixel wide curve, are: Precision, skeleton-based Recall and skeletonbased balanced Accuracy.
The Precision is calculated as follows: where TP (true positives) represents the number of pixels belonging to model M within the MTA delineation and FP (false positives) represents the number of pixels belonging to model M outside the MTA delineation. Following this idea, a skeleton-based Recall can also be defined, as shown in (15), to indicate the ratio of correct positive pixels out of the total positive pixels that a perfect model should contain.
where P represents the number of pixels in the skeleton of the MTA delineation, which will be taken as the ideal number of pixels that a one-pixel wide model should contain. Finally, a skeleton-based balanced Accuracy is considered to measure the general performance of the model, using the following definition: with TPR = TP/P, and using N (skeleton negatives) as the number of pixels that do not belong to the skeleton of the MTA delineation.

Implementation Details
The proposed method was evaluated using the dataset DRIVE with MTA manual delineations.
Firstly, for the segmentation part, The U-Net network was trained with 200,000 random patches with size 48 × 48 pixels, from which 180,000 were used for training and 20,000 were used for validation. The optimization process was performed for 150 epochs, using the Stochastic Gradient Descent (SGD) optimizer with minibatch size of 32 patches.
Secondly, for the spline curve construction in the numerical modeling part, a configuration of five points with quadratic functions has been chosen.
Finally, since the numerical modeling part contains a stochastic component for point selection and model construction, 30 independent executions were carried out and the averaged results are reported.
All the experiments were executed using a computer with an Intel Core i5 2.4GHz processor and 12GB of RAM, excepting the segmentation step, for which a Nvidia Tesla K80 GPU with 12GB of RAM was used. The source code is freely available (https://github. com/dora-alvarado/robust-MTA-detection-modeling, accessed on 11 April 2022).

Comparative Analysis
The proposed method was evaluated using the test set of the dataset DRIVE, consisting of 20 images with size 565 × 584 pixels. In Figure 5, some examples of the MTA numerical modeling are presented. The performance analysis was carried out considering four Hough-based state-of-theart MTA detection methods: the Gabor-based enhancement combined with a Hough detector (Gabor+Hough) [7], the General Hough detector (General Hough) [48], the parabola detection algorithm from MIPAV software (MIPAV) [49], and the hybrid UMDA+SA parabola detector (UMDA+SA) [11].
A comparison in terms of closeness between the numerical model and the MTA manual delineation is shown in Table 1. The proposed method reaches the best value in MDCP and Hausfford Distance, more than three pixels and six pixels away from the second best, Gabor+Hough [7], respectively. A comparison for the MTA detection task has also been made using Precision, skeletonbased Recall and skeleton-based balanced Accuracy. As shown in Table 2, the proposed method obtains the best values in the three metrics, doubling the performance of the second best (UMDA+SA [11]) in Precision.
Through a qualitative comparison presented in Figure 6, it can be inferred that the difference in performance lies in the ability of the methods to adjust to a non-symmetric MTA shape. The proposed method shows a robust behavior in these scenarios, while the rest of the methods diverge from the manual delineation due to its parabolic foundation.  Finally, a comparison considering execution time is reported in Table 3. The proposed method obtains the second best value, only surpassed by the hybrid method UMDA+SA [11], and with an execution time 20 times faster than the third best, Ga-bor+Hough [7]. Table 3. Execution time comparison of MTA detection and modeling using the test set.

Conclusions and Future Work
This paper proposes a new method for automatic MTA detection and modeling. The segmentation step contributes to the computational efficiency of the proposed method, reducing the search space to only blood-vessel-related pixels. Unlike the previous works, based on semi-elliptical parabolas, the proposed method relies on a piecewise parametric function with the ability to adequately represent both symmetric and asymmetric MTAs. Through a weighted-RANSAC scheme that takes advantage of a priori knowledge about the MTA characteristics, the algorithm makes a robust selection of points to build the model. The inclusion of blood-vessel width and foreground-location estimations for inlier counts promotes selecting a model built with high probability MTA points. The method has proven to be robust and efficient in the MTA modeling task, obtaining a Balanced accuracy of 0.7067, MDCP of 7.40 pixels, Hausdorff distance of 27.96 pixels, and an average execution time of 9.93 s per image. These numerical results have also shown that the proposed method is suitable for implementation in systems that perform computer-aided diagnosis in ophthalmology.
A future direction of this work may be to determine the MTA openness for the presented piecewise-parametric model with the aim of classifying ophthalmological alterations. The first approach to this task would be to quantify the area-under-the-curve and slope variations for the function.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: