Log In Sign Up

Fully Automated Tree Topology Estimation and Artery-Vein Classification

We present a fully automatic technique for extracting the retinal vascular topology, i.e., how the different vessels are connected to each other, given a single color fundus image. Determining this connectivity is very challenging because vessels cross each other in a 2D image, obscuring their true paths. We validated the usefulness of our extraction method by using it to achieve state-of-the-art results in retinal artery-vein classification. Our proposed approach works as follows. We first segment the retinal vessels using our previously developed state-of-the-art segmentation method. Then, we estimate an initial graph from the extracted vessels and assign the most likely blood flow to each edge. We then use a handful of high-level operations (HLOs) to fix errors in the graph. These HLOs include detaching neighboring nodes, shifting the endpoints of an edge, and reversing the estimated blood flow direction for a branch. We use a novel cost function to find the optimal set of HLO operations for a given graph. Finally, we show that our extracted vascular structure is correct by propagating artery/vein labels along the branches. As our experiments show, our topology-based artery-vein labeling achieved state-of-the-art results on multiple datasets. We also performed several ablation studies to verify the importance of the different components of our proposed method.


page 3

page 5

page 7

page 9

page 10


Fully Automated Artery-Vein ratio and vascular tortuosity measurement in retinal fundus images

Accurate measurements of abnormalities like Artery-Vein ratio and tortuo...

Joint segmentation and classification of retinal arteries/veins from fundus images

Objective Automatic artery/vein (A/V) segmentation from fundus images is...

Simultaneous segmentation and classification of the retinal arteries and veins from color fundus images

The study of the retinal vasculature is a fundamental stage in the scree...

Efficient Kernel based Matched Filter Approach for Segmentation of Retinal Blood Vessels

Retinal blood vessels structure contains information about diseases like...

Iterative Deep Learning for Network Topology Extraction

This paper tackles the task of estimating the topology of filamentary ne...

Curvature Integration in a 5D Kernel for Extracting Vessel Connections in Retinal Images

Tree-like structures such as retinal images are widely studied in comput...

1 Introduction

Retinal fundus images allow ophthalmologists to diagnose a variety of ocular and cardiovascular diseases, including diabetic retinopathy (DR) [retinal_imaging_and_analysis-5660089], glaucoma [economic_analysis_pmid10868871], age-related macular degeneration (AMD) [AMD_severity_scale-Ferris2005], and the likelihood of stroke [retinal_signs_and_stroke_baker_retinal_2008, Kipli2018]. These diagnoses are based on tell-tale features in the fundus image that are correlated with a higher likelihood of certain diseases. For example, a higher cup-to-disk ratio in the optic nerve is correlated with a higher likelihood of glaucoma [cup-to-disc-glaucoma-BOCK2010471]. Similarly, the artery-vein ratio can be used to predict a patient’s risk for diabetic retinopathy and hypertension [retinal_imaging_and_analysis-5660089, hyper_pakter2001detection, av_ratio_10.1016/j.amjhyper.2004.10.011, retina_vascular_calliber_ikram_2013]. Other diagnostically useful features include vessel tortuosity, bifurcation and branching angles and the presence of exudates [10.1093/eurheartj/ehm221, 10.1167/iovs.03-1390, doi:10.1117/12.708469, Lue31, retinal_imaging_and_analysis-5660089, exutates_DR-JOSHI20181454, turtuosity_use_HART1999239].

Figure 1: Flowchart of the whole pipeline: We start from a fundus image then we generate Vessel Probability map(Vessel prior, fig. 2.c), Artery/Vein Probability Map(AV Prior, fig. 2.e) from two separately trained U-Net CNNs, we use multilevel skeletoniziation(algo. 1) to produce union graph of the vasculature(fig. 3 a & b) which is further pruned significantly with a unique graph contraction technique(algo. 2, fig. 3 c & d). Them we use parallel Dijkstra’s shortest path algorithm to extract a un-directed topology graph, which is assigned direction by our Flow assignment algorithm to obtain a directed graph() as in fig. 3.e. We then map the labels to as shown in fig. 6.a. Then we perform series of High Level Graph operation(HLO) to minimize topology cost function Eq. 3. We then perform a simple A/V label propagation algorithm to show that our topology graph is correct using different measures explained in experiment section. Ablation: We perform heavy ablation studies(shown by blocks in light blue color) by replacing vessel prior with manual vessel segmentation, and AV prior by AV manual segmentation, alternating one another such that at a given point we have at least one prior in the pipeline. The pipeline works best when we have good priors(Vessel segmentation and AV label priors), as such we can always improve our priors and come back to obtain better performance.


Medical image analysis tasks mostly fall in classification, and segmentation tasks as we need to identify some region with some abnormality to detect diseases. As such, we use Fully Connected Convolutional Neural Network (FCNN) to perform segmentation task. We assign a class label (0, 1, …,

) to each pixels in segmentation task. U-Net, a multi-level encoder-decoder based architecture is uses as a backbone for major medical image segmentation task [Ronneberger2015UNetCN]

. It has a bottleneck part on each level that helps ignore noise and focus on information that matters, also allowing the skip connection to keep any information lost in the process. It uses sliding window approach to make such architecture to work with any size images so we have used U-Net as the backbone for our intermediate feature extraction tasks as in figure


, where we use two U-Nets for vessel probability-map assignment and Artery/vein probability-map assignment. Usually, we use a threshold to binarize such probability-maps to obtain a vessel mask. Using a fixed mask can lead to a very annoying trade-off problem of precision and recall: very low threshold captures too much noise, whereas a high threshold misses a lot of fainter vessels. More importantly, only a binary mask is not enough to measure features like tortuosity, branching factors, vessel crossing anomaly, Artery/Vein ratio and more. Thus in this paper we aim to extract the entire vasculature as a graph topology such that we can traverse to any part of the vascular structure and take all sorts of before-mentioned measurements, also run complex structural analysis on it

[graph_concepts_PATTON200699, retina_vascular_calliber_ikram_2013].

In order to capture the complete vasculature, we use a novel idea of multilevel-skeletonization (3.1) where we use a set of thresholds ranging from low-to high and capture the complete vascular structure. Then, we convert the combined skeleton generated from binary images of each threshold and convert that into a lattice graph(1) called union-graph. This union graph can be seen in figure 3.a,b where the cyan colored nodes are branch-nodes(), golden nodes are called end-nodes(). Naturally the branch nodes should only appear in the branching locations, however, we can see the union-graph being dense. Thus, we use a novel pruning technique called union-graph-contraction (2) to get rid of significant amount of nodes and edges keeping the overall vascular structure intact( 3.c). Then we use 2-step Dijkstra’s shortest path algorithm(4, 4) to extract the clean graph structure of the vasculature as shown in figure 3.e. We then assign direction7 on this contracted graph along with the concept of Pseudo Optic Nerve Head Assignment to detect disconnected regions as shown in figure 3.f. The extracted directed topology graph opens up numerous potential research paths for downstream analysis and diagnosis tasks. For a instance, in this paper, we used the topology to propagate Artery/Vein label to obtain better A/V segmentation mask which is used to calculate a crucial bio-marker called Artery/Vein ratio that reveals arteriolar narrowing liked to hypertensive retinopathy, stroke and coronary artery disease [av_ratio_10.1016/j.amjhyper.2004.10.011, retinal_signs_and_stroke_baker_retinal_2008]. Similarly, one can easily run analysis like branching factor analysis, turtuosity, crossing anomaly and so on as shown in figure 1.

We used two priors to do A/V label propagation;A/V probability obtained from U-Net, directed graph topology obtained by flow-assignment algorithm(7). There are some inherent limitations stemmed from the skeletonization algorithm in union skeleton algorithm (1) that shifts some crossing either upward or downward introducing unnecessary sink and source nodes as shown in white-coloured nodes in figure 9 a & b respectively. Section 3.2.3 list a handful of simple high level operations(HLO) that we can preform on the directed topology graph to fix such problems. We simply minimize a vessel topology cost function 3 by performing series of such operation as mentioned in figure 7, 8. Once we have a low cost topology graph, we simply propagate Artery/Vein labels by using simple propagation strategy discussed in section 3.3. We also did very crucial ablation studies that reflects on how the quality of priors (Vessel probability and Artery/Vein) affects the overall quality of the topology, and how it affects the downstream tasks like A/V label propagation.

Figure 2: a). Original images b). Manual vessel segmentation c). Vessel prior obtained from U-Net d). AV prior generated from U-Net


2 Related Work

Fundus images widely used for vesicular and eye related disease detection. A good vascular topology can lead to more accurate disease detection. Over the years many researches tried different approaches to tackle this challenging problem. Some researchers applied graph based approaches to solve this problems [Amil2019], [8754802], [7120990]

. In recent years many researchers applied deep learning technique for both segmentation

[6848752], [HEMELINGS2019101636], [GIRARD201996] and classifications [XU20173], [8999616], [Hu2021] of fundus images. In the following paragraphs we review some of these work in graph topology, segmentation and classification of the fundus images.

2.1 Graph topology

In [Amil2019], Amil et al. applied the structure of the blood vessels to detect patient with diabetic retinopathy and glaucoma. First, use automatic unsupervised segmentation algorithm to extract the tree-like graph. Then quantify structural differences between the graph of healthy and non-healthy patients. Fractal analysis used to characterize the extracted graphs. The paper claimed there is a significant differences between healthy and mom-healthy images. Zhao et al. in [8754802] used a dominant set clustering and formulated the retinal blood vessel topology estimation and A/V classification in the pairwise clustering system. Image segmentation, skeletonization and significant nodes identification used for graph construction. Also, the invert Euclidean distance applied between two end points. Classification into A/V applied based on intensity and morphology after constructing the vascular network. In the previous work Estrada et al. in [7120990] applied graph-theoretic framework to differentiate arteries from veins. The use of underlying vessel topology made a better classification for small and midsize vessels. Applying the global likelyhood model explore the possible solutions relevant to project vessels. The model is capable of analyzing the entire vasculature. The current proposed model is an extension and fully automated version of this prior work.

2.2 Segmentation

Roychowdhury et al. in [6848752]

applied three-stage algorithm for blood vessel segmentation. In the first stage a binary image extracted from the high-pass filtering on the green plane of fundus image. Also another binary image of vessel region extracted from the morphologically reconstructed enhanced image. If the region is common in both binary images, it will be extracted as a major vessels. In the second step, by using a Gaussian mixture model (GMM) the remaining pixels in two binary images are classified. The classifier use a set of eight features which are extracted based on the pixel neighborhood, first and second order gradient images. The major portions of the blood vessels and the classified vessel pixels are combined in the third step. Roychowdhury et al. claim less dependent on training data, less segmentation time and consistent vessel segmentation accuracy on normal images in their algorithm. Hemelings et al. in

[HEMELINGS2019101636] applied deep learning architecture for A/V segmentation. To be more specific U-net architecture applied for segmentation. The model tested on DRIVE and HRF datasets.Girard et al. combined deep learning and graph propagation in [GIRARD201996]

for A/V segmentation. The convolutional neural network (CNN) used to jointly segment and classify vessels into arteries and veins. After the initial CNN labeling, a graph theory used to calculate the cost of linking the pairs of branches which represent vessels. A minimum spanning tree use to simplify the graph and improve the efficiency.

2.3 Classification

Xu et al. used automated computer-aided method to predict A/V classification. In order to reduce the differences in feature space, the Intra-image regularization and inter-subject normalization applied. This method tested on DRIVE Dataset[XU20173].An attention guided U-net with atrous convolution(AA-UNet) applied by [8999616] to segment the vessel and non-vessel pixels. They use AA-UNet to generate the attention mask and applied as a weighting function to increase the network attention toward the vessel region. Atrous convolution instead of ordinary convolution in feature layer to increase the receptive field and decrease computation. In [Hu2021], Hu et al. applied vessel-constraint network(VC-Net) to enhance A/V classification. This method utilizes the information of vessel distribution and edge. VC-Net uses a vessle-constraint (VC) module which combines local and global vessel information. The model generates a weight map to constrain the A/V features which enhance the edge and end feature of blood vessels. It also suppresses the background-prone features. In order to to extract the blood vessels information and improve the feature extraction capability and robustness of the model, a multiscale feature (MSF) module used.The VC-Net get vessel segmentation result simultaneously.

3 Methodology

In this section, we describe the three major steps of our fully automated vessel extraction method: (1) topology extraction with multilevel skeletonization, (2) topology estimation, and (3) artery-vein label propagation. Table 5 lists the symbols used throughout the rest of this paper. Super- and subscripts are only used when necessary to differentiate between different stages of the algorithm.

3.1 Topology Extraction with Multilevel Skeletonization

The first step is to extract vascular topology and measure different important factors like (Artery to Vein )AV-ratio, turtosity, bifurcation mentioned in the section 1. Its crucial to correctly segment vasculature in order to compute usable such measures in the original image, thus our method works on the original image size without having to resize. Funduscope, an imaging device to capture fundus images comes in different size and quality, can have different resolution images and we have shown that our technique works in variety of datasets end-to-end with no to minimum parameter optimization. Let (fig. 2 c) be the likelihood map of a fundus image (fig. 2 a) generated by a trained neural network.

Figure 3: a) Undirected dense graph representing vessel path (algo. 1). b). Zoomed dense graph to show nodes and edges c). Pruned graph with significantly less nodes and edges but with intact vasculature (algo. 2) d). Zoomed pruned graph e).Result of Parallel Dijkstra shortest path applied on pruned graph (algo. 3 and 4). f). Direction assigned to graph obtained from Dijkstra’s shortest path algorithm with flow assignment algorithm(algo. 7). For better visualization, we have only shown the median node of each segments and branch nodes. We only use vessel likelihood to generate such directed graph, however for High Level graph Operations(HLO) and label propagation we map the Artery/Vein prior to this final directed graph and go from there.


3.1.1 Extraction of Likelihood-Map

We use Convolutional Neural Network on the green channel to get vessel likelihood map as . Specifically we use a U-Net architecture [Ronneberger2015UNetCN] with stochastically weighted loss [deeddyn_10.3389/fcomp.2020.00035] to train k-fold for DRIVE [drive_dataset] and WIDE [wide_dataset:6987362]. The stochasticity helps detect finer vessels with ease[deeddyn_10.3389/fcomp.2020.00035].

3.1.2 Extraction of AV-Prior

We extract likelihood for Artery/Vein (noted ) for each pixel by k-fold cross validation with U-Net. Such pixel-wise classification is inconsistent to actually use for measures like AV ratio, bifurcation because of too many mis-classification and noise as seen in figure 2 e. Thus, we use it as a prior to topology optimization and AV classification.

3.1.3 Multilevel Skeletonization

We then use a range of thresholds given by on , which ranges from with step size of to extract multilevel-skeleton(algo. 1); This skeleton is converted to a dense lattice graph called union graph (3.a). We further used a novel pruning technique(algo. 2) to further contract the graph nodes/edges to significantly reduce the size (fig. 3.c).

= Vessel likelihood map, where each pixel range from .
= 0 initialized array being same size as the image.
for t in range(255, , -p) do
       = ;
       = ;
       = );
end for
Algorithm 1 Multilevel skeletonization of : function sets every pixel in to 0 if less that t, else 255, can be any skeletonization morphological operation, and is pixel-wise maximum. Now, union graph is simply a lattice graph of , where each node is a pixel with value 255 in , and two nodes have an edge if they are 8-neighboring pixel.

Each node of graph has three attributes; Likelihood of corresponding pixel , Color of pixel in Lab space of original fundus image, and Shortest distance to background in binary image with lowest threshold in , . We further reduce number of nodes in union graph() by using the following contraction technique in algorithm 2. It fully preserves the vasculature however reduces the graph size(nodes and edges) significantly as shown in fig. 3.c.

= Union graph.
= .copy()
= Branch nodes in sorted in increasing order of vessel likelihood.
while  is not empty do
       = pop();
       if not in : continue;
       for n in -  do
             if not in : continue;
             if : ; ;
             if in : ;
       end for
end while
Algorithm 2 Union graph contraction to get the lightweight contracted union graph

3.1.4 Dijkstra shortest path algorithm for complete vessel path tracing

We use Parallel Dijkstra Shortest path algorithm to estimate vessel path on the graph

. Being lightweight and parallel, the running time is greatly reduced. We use the geometric mean of the node attributes mentioned in section

3.1.3 to weight the edges as below:


, where

  • = 255 - , where is the vessel likelihood of the pixel .

  • = color similarity metrics (deltaCie2000(n1, n2) in Lab color space).

  • = +

We use two pass of this algorithm (algo. 3 and 4) because shortest path algorithm rely on greedy approach to be cost effectiveness and could always favour some paths over the other. This introduces un-visited vessels, leading to disconnected paths. In order to fix such issues we, first detect such lumps of un-visited regions, normalize edge weight based on distance to the current extracted topology(alg. 5) and rerun the algorithm in such regions only (4).

= ;
for c in sorted desc do
       = ;
       src = ;
       if not in : src = [0];
       for parallel t in  do
       end for
      ; ;
       if : ;
       for parallel t in  do
       end for
end for
Algorithm 3 Parallel Dijkstra’s’ shortest path algorithm for vessel tracing: This algorithm, with some fine tuning with algo. 4, takes in the contracted graph from algo. 2, fig. 3.c and outputs a undirected topology graph as in fig. 3.e
= ;
= weighted contracted union graph;
for i in  do
       = ;
       for  in  do
             ; ;
             = ;
             ; ;
             if : continue;
             = set();
             for p in  do
                   if in : ;
             end for
            (, ;
             for seg in  do
             end for
            for b in  do
             end for
            if : ;
       end for
end for
Algorithm 4 Fine-tune Dijkstra’s algo. 3: The undirected graph extracted by shortest path topology extraction algorithm is not vessel topology aware, which means it follows the shortest path on a given cost, but some images can vary in brightness, contrast, some will have medical abnormalities leading to always follow a cost-effective path creating disconnected vessels. In order to fix that, this algorithm readjust weights() in such regions where a significant portion of the contracted graph is not visited by the algo. 3. It is parameterized by as in algo. 5
       = - ;
       for n in  do
             = ;
       end for
      for , in  do
       end for
Algorithm 5 - Re-initialize edge weights based on how far the left nodes are form the currently extracted topology: The parameter controls how the edge weights are influenced by the distance of two edge nodes from the current topology.

3.1.5 Graph cleaning

We clean obvious noise in the graph as folows:

  • We use a radius of to contract optic disk as in figure 3.f where we ignore nodes withing the disk.

  • Remove self loop edges and isolated nodes.

  • Remove leaf segments with less that nodes.

  • Remove small cycles with nodes.

  • Disconnect corner connections(abrupt turns).

  • Node fill for uniform nodes distribution.

  • Smooth graph path to center on vessels.

3.2 Topology Estimation

This step converts the topology obtained from previous steps to a directed graph with approximated blood flow such that we can optimize the topology. Without direction assigned the estimating correct topology is NP-hard, thus our flow assignment algorithm 7 is robust to assign approximate blood flow.

Figure 4: Pseudo ONH obtained by our flow-assignment algorithm(algo. 7) The cyan node being the true optic nerve head (ONH), whereas yellow are the end nodes(degree = 1). The size of the node is the proportion of number of flow ends visits to that node in flow-assignement algorithm. In short, we traverse from the end points and branch node neighbors until we hit either ONH or another end node. We canW see that most of the visits are accumulated at true ONH, also in nodes which are the disconnected points in the vessels.


3.2.1 Flow assignment

This algorithm takes an extracted undirected graph by Dijsktra’s algo. 3 & 4 and assigns direction to each edges as in figure 3.f. It is a crucial step as direction is one of the most important prior for proper topology extraction along with Artery/Vein prior labels. It is based on an elegant recursive tracking algorithm that uses and branching trait of vasculature (algo. 7). Branching from a parent vessel tend to branch with and branch away from main vessel coming from optic nerve head(ONH). It also has a momentum factor involved where collective weak flow converges to a stronger one just like river networks. The core of this algorithm is the route function (alg. 6). It traces vessels from end nodes() and branch nodes() and count number of visits on each checkpoint nodes(–end nodes and optic nerve head(ONH)) of an undirected graph . The route function takes a branch-forward method that determines what metrics to use to choose the path in a branch. One can choose to use any branch-forward measure, however we choose the following:


where is the straightness of segment when forwarding to another segment of a branch s.t. (fig. 5).

= any un-directed graph.
= any branch node in .
= any adjacent node of .
path = route{
       if in : return ;
       = ;
       return ;
Algorithm 6 Function Path route. Where, the branch forward function is given by eq. 2. It uses a vessel straightness measure(fig. 5) based on vessel orientation. Note that we can replace the branch forward method with any desired method.

We use the branch-forward method shown in eq. 2 for assigning direction with the flow assignment algorithm. This method relies on a regularized version of straightness measure based on angle between two segments parametarized by , and such that small mis-alignment of neighboring nodes does not lead to wrong choice and yields strong confidences(fig. 5).

We can see in figure 4 where node size corresponds to number of visits on each checkpoints–most of the convergence is on optic nerve head (ONH). More importantly, also in the places where a vessel is disconnected due to some pathological condition.

Figure 5: Default branch-forward method on each branch nodes based on multi-part orientation of branch segments(). The checkpoint of interest here is and we aim to forward this flow from segment through branch .
= .
= {} .
for  in +  do
       for  in  do
             = ;
       end for
      for  in  do
             = ;
             = ;
             = + ;
             = + ;
             if : = ;
       end for
end for
Algorithm 7 Flow assignment and pseudo-optic nerve head: This algorithm relies on the extracted undirected graph as in fig. 3.e and generate directed graph 3.f. We traverse towards from each end nodes(, ), and from neighbors of every branch node (), where is the undirected graph obtained from 5. We then traverse with the routing algorithm 6until we hit a checkpoint node (). We count number of hits on such nodes depicted by node size in fig. 4.
Figure 6: Prior Mapped to flow-graph obtained from algo. 7. i) Flow obtained from vessel prior(from U-Net), and AV-Prior(U-Net). ii) Flow obtained from vessel-Prior and AV-Ground Truth. a) End nodes(golden yellow, ), branch nodes(cyan, . b) Source branch nodes(white), c) Sink branch nodes(white), d) Graph after High Level Graph Operations(HLO) + A/V label propagation. e) Artery/Vein segmentation mask generated from propagated topology graph by mapping it to a binary segmentation and using nearest neighbour label propagation.


3.2.2 Tree topology cost function.

We use an intuitive cost function to optimize the topology:


, where is the un-directed graph of , is crossings cost that looks if illegal crossings are made, restricst the correct distribution of A/V in branches, is branchness cost that respects the true nature of vessel branching, and cost keeps the blood flow between A-A/V-V in check.


, where is number of arteries segments connected to branch , is the degree of branch , is the normalization factor enforcing number of arteries and veins in a crossing should be same when degree of that crossing is 4. And should try to be in pairs when more. We use (instead of 2) to make sure that we do not penalize natural crossing for cases with .


, where is the average artery expectation of nodes in neighboring segments of branch , and likewise for veins. It enforces that on each branch we make sure the distribution of A/V is near discrete(either 0 or 1). Any violation to this triggers both and .

Similarly, algo. 8 is the cost incurred when any violation to branching is made–such as incoming vessel in a crossing must have an outgoing vessel, or out-degree should always be more than in-degree. Also, when there are multiple vessels crossing, the segment that share a minimum branch forward cost () as in fig. 5 should be preferred.

       if : ;
       if : ;
       return ;
Algorithm 8 Branchness cost for each branch in a directed graph , where is the in/out degree at branch of the corresponding directed graph .

The cost is incurred when the branch-forward done by (algo. 8) does not align with a strong prior, which means the two favourable segments now has very weak compatibility(eg. one is close to Artery, other to vein). For that we introduced a branch-forward method that has the A/V prior factor added along with the orientation(straightness) as mentioned in fig. 5 as follows:


, where s.t. . For all in-neighbors of a branch, a high cost is incurred if the most straight outgoing has different label.

3.2.3 High Level Graph Operations (HGO) to enable blood flow

Graph topology estimation is a NP-hard problem [wide_dataset:6987362], however a good directed graph from section 3.1, and artery-vein prior from UNet from section 3.1.1 makes the optimization problem much easier-we are using the prior mapped to graph and minimizing the equation 3. In order to do so, we have defined a finite set of high level graph operations (fig 7, 8). It is important to understand how such scenarios occurs so it demands such HLO in the first place. When two vessel(Artery and Vein) cross each other, most path extraction algorithm including ours tend to shift the crossing either up or down creating two 3-degree source and sink nodes instead of actual 4-degree branch node(fig. 8 HLO1.i.), and the defined HLO operations are a reverse of such malformations to bring back the crossing to its original state.

We use the directed topology graph mapped with A/V labels from to identify fault points; source and sink. On these nodes, we perform HLO and maintain a heap of edited graph with the graph cost (eq. 3). The top of the heap is considered the next best topology state only if the cost difference with the current best state is more than . We want the optimization algorithm to avoid any noise and only keep robust and correct operations in the optimization path. If any true operations are missed because of weak prior, the algorithm being iterative, will eventually propagate to make the necessary jump. This cost is subject to only local changes thus won’t detect any misconnection far off from the source/sink of operation, so each time after a HLO, we propagate labels from the shift destination nodes(fig. 7 HLO1.i, HLO2.i, HLO3.ii, HLO4.ii, and HLO5.i) before calculation graph cost and push to the heap. This makes sure that our operation is a valid one locally and in the global vasculature context.

Source nodes: These are the branch nodes that are bifurcated in the directed graph (fig. 8, 9.a). We use the route method (algo. 6) on un-directed graph of at two outgoing vessel with the flavoured forward algorithm (eq. 6) instead of default (fig. 5). On such bifurcation points which gives us two vessel paths & . A true source is when two vessels are different from one another(i.e. one artery, other vein) because arteries and veins should not cross each other, as such the can never bifurcate from that point. To avoid propagation that violates this fact, we calculate a propagation score (eq. 8) and select nodes only greater than a threshold prop_score_thr. In addition to this threshold based selection, we consider a bifurcated branch node source if the two outgoing sub-graph intersected down the line somewhere. This is a necessary constraint because as mentioned before arteries and veins never cross each other thus if the two paths meet, both of them has to be different labels. Even if it is because of wrong prior, it is actually better to mark this node as a source because our optimization iteration labels more confident cases first, thus improving the prior to make better decisions.


, where , and , where is just a scaling factor for adjustment if necessary.

Figure 7: High Level Operation Shift up and Shift Down from sink.


Figure 8: High Level Operation Shift Down from source.


Figure 9: Source and sink nodes in actual fundus image patches


Sink nodes: Sink nodes are branch points in directed graph where two vessels converge to a node(fig. 7, 9.b). Similar to source nodes, a true sink node is actually the one where two incoming vessels are different. We also use propagation score based threshold to select sink nodes as it makes sure we work on confident propagation first, improving expectation on rest of ambiguous decisions.

High Level Operations: There are four basic operations on sink nodes as in figure 7: Shift-Up, Shift-Down, Reverse-Shift-Up, and Reverse-Shift-Down. One can notice that we have to make a choice for Reverse operation in fig. 7 that which segment to reverse because there are two choices(two inward vessels in the sink node, fig. 7.3,4 i). A close inspection reveals that we can just pick the segment with highest cost(eq. 4). There are two operation in source nodes; one of them is Shift-Source-Up as shown in figure 8, and there could be another one similar but Shift-Source-Up which is found to be very rare in our dataset so we discarded in our experiments but one may include with minimal change if necessary.

High Level Graph Operation on sink: We perform the before-mentioned four basic operations on each sink nodes sorted in ascending order by the distance to ONH and push to a min heap–The top of the heap is the topology with minimum cost. Working with sink nodes that are closest to ONH first makes more sense because routing from them lead to ONH through other sink nodes that are already worked with, thus minimizing the chance of worst case operations. We can see sink nodes in figure 6.c, whereas fix applied with HLO operation in 6.d. One can notice that a sink node must be accompanied by a source node in order to make a correct shift UP/DWN operation.

HLO1:ShUP We simply merge the outgoing segment from the sink down to the middle node of that segment. This in turn squeezes the outgoing segment and its end branches to a single branch as in fig 7.i. We intuitively called this shift-up operation because we are shifting the sink upwards in the direction of the flow. Also, we start working from the sink that is near to the optic nerve head.

HLO1:ShDWN Shift down operation is similar to shift up operation but in reverse direction fig. 7

Reverse The flow assignment algorithm (alg. 7) depends on how good the vessel prior is. So there could be cases where a segment is assigned wrong direction, with only a reverse operation it would be fixed like in fig. 7 HLO3.i, HLO4.i. It should be noted that not all HLO operation lead to a good optima. The core of such operation is fueled by our graph cost measure in eq. 3.

High Level Graph Operation on source: As mentioned before, source and sink nodes come together in majority of the cases, however there are some cases where source node appears by themselves–Usually such locations forms where a vessel passes through the branching of another vessel like in fig. 8 HLO5. In theory, we can find both shift up and shift down operations in such locations, but our closes inspection in dataset reveals a handful cases of shift down, and almost negligible cases of shift up such operations, so we only have listed HLO5:ShDWN-SRC (fig. 8). Unlike sink node, source nodes are always worked with from the ones furthest from the ONH for the similar reason as sink nodes.

3.3 AV Label Propagation

We have explained the vasculature extraction process in section 3.1. In order to validate the claim that topology extracted from our technique is correct, we begin with av-pmap and the directed graph (, and for corresponding undirected graph) extracted from section 3.1, and use a simple iterative propagation algorithm to propagate A/V labels correctly. As such, our results have consistently shown improvement in A/V scores as shown in Experiments and Results section (Sec. 4). One of the challenges of graph based label propagation is that a minute false-propagation could affect all the other regions depending on it. For example, a false main vessel label could end up labeling all the branch vessel of it falsely. Thus we have a concept of stop-points: These are nodes from where we avoid any propagation further such as a crossing but where one of the neighboring segment is undetectable because of some condition, and this 4-neighbor branch(crossing) end up being a 3-neighbor branch(Source node) and one of the neighbor being different label. Another example is, when to directed paths from source nodes meet down the tree violating A/V crossing constraint. We can enforce to avoid routing beyond such stop-points by just adding these stop-nodes to the checkpoint nodes of the corresponding graphs() in the routing algorithm(algo. 6). We use two propagation points:

Propagate from sink We start propagating from the sink nodes that are closest to the optic nerve head (ONH). This minimizes the probability of false propagation because all the paths we are about to propagate labels to are already worked with. We simply set the probability of each nodes of the two incoming path obtained by route algo. 6 to the corresponding average expectation. Once we are done with all sink nodes, we do it iteratively until the graph cost eq. 3 improves no more.

Propagate from end nodes. After propagating from sink, we propagate from end nodes to the ONH. We use the same routing method algo. 6 from each end nodes and replace each path by the average expectation of the vessel. The average expectation works well in this scenario because A/V priors obtained by neural network are near discrete(close to 0 or 1). We also do this with multiple passes over all end nodes and stop when no more desired improvement is noticed in the graph cost.

4 Experiments and Results

We have done experiments in three datasets DRIVE [drive_dataset], WIDE[wide_dataset:6987362], and INSPIRE[inspire_dataset-tang2011robust]. This work uses two different priors: vessel likelihood() to extract topology in which each pixel is assigned a probability of being vessel, and background. A/V likelihood to propagate AV-Label() in which each pixel is assigned 3 probabilities as a distribution—artery, vein, and background. in subscript specify likelihood prior obtained by a trained U-Net, whereas represents ground truth as a prior, whereas in superscript represents the state before A/V propagation, and (also bold) represent score after A/V propagation. We have used ground truth of either or to ablate the two priors to see the ablation effect of quality of priors, however both ground-truth never come together in one experiment.

4.1 Transfer learning for INSPIRE dataset

For INSPIRE dataset[inspire_dataset-tang2011robust], we do not have a vessel segmentation, or AV segmentation mask, however we do have a labeled topology with sparse nodes specifying the topology, lets call these set of nodes(pixels)

. We trained a model using U-Net architecture on other datasets(DRIVE, HRF, CHASEDB) by resizing images and ground truth to approximately match with that of INSPIRE dataset, also normalizing with the training mean and standard deviation. Then we generated vessel likelihood map for all images in INSPIRE dataset. Thus in our

combination column in tables 3, 4 we do not have an entry for . We then binarize the likelihood map, and for each vessel pixel in this binary image we assign AV label to whichever node in is closest(Nearest Neighbor). We treat this AV labeling as AV-segmentation ground-truth . We then train U-Net with cross validation to generate AV likelihood map for all images.

4.2 Topology Extraction by AV propagation: ablation study

We have calculated scores with the ground truth in three different ways(listed in different columns groups with different superscripts—, , and ):

4.2.1 Node-node comparison

is the score of nodes of graph from section 3.1 labeled based on corresponding prior in combination and mapped with ground-truth. As such, we label each node of with either Artery or Vein. As mentioned before, the combination in bold is after we apply our AV propagation algorithm.

4.2.2 Common pixels comparison

We also asses by AV labeling a binarized vessel likelihood map by using Nearest Neighbor approach as mentioned in section 4.1. We the take only the common pixels between such labels and ground truth, and calculate score. Here, the pixel under consideration will only have either Artery or Vein label since taking intersection omits background.

4.2.3 Complete AV-label comparison

In this approach, instead of only considering common pixels, like section 4.2.2, we compare complete pixels in the image-we will have 3 possible labels for each pixel—artery, vein, and background. Then we use the manual AV ground-truth to calculate score. We keep record of this score before and after applying propagation algorithm.

Precision Recall F


0.9181 0.9898 0.9708 0.8872 0.9895 0.9708 0.8941 0.9896 0.9708

0.9325 0.9638 0.9684 0.9274 0.9643 0.9684 0.9287 0.9638 0.9684

0.7814 0.8738 0.9649 0.7822 0.8723 0.9649 0.7784 0.8710 0.9649

0.8171 0.8806 0.9653 0.8180 0.8798 0.9653 0.8143 0.8781 0.9653

0.7780 0.8815 0.9625 0.7808 0.8795 0.9625 0.7754 0.8783 0.9625

0.8140 0.8888 0.9631 0.8174 0.8877 0.9631 0.8126 0.8863 0.9631

Table 1: Topology estimation results on DRIVE dataset:

Precision Recall F


0.9117 0.9947 0.9753 0.8816 0.9947 0.9753 0.8858 0.9947 0.9753

0.9319 0.9670 0.9726 0.9299 0.9663 0.9726 0.9302 0.9665 0.9726

0.8122 0.8850 0.9721 0.8114 0.8825 0.9721 0.8102 0.8822 0.9721

0.8739 0.9127 0.9736 0.8724 0.9111 0.9736 0.8715 0.9107 09736.

0.7982 0.8955 0.9681 0.7989 0.8931 0.9681 0.7967 0.8927 0.9681

0.8305 0.9042 0.9685 0.8311 0.9031 0.9681 0.8285 0.9026 0.9685

Table 2: Topology estimation results on WIDE dataset:

Precision Recall F


0.9412 0.9873 0.9921 0.9369 0.9869 0.9921 0.9370 0.9871 0.9921
SEG + AV 0.9367 0.9447 0.9890 0.9353 0.9445 0.9890 0.9350 0.9439 0.9890

0.8105 0.8696 0.9839 0.8092 0.8695 0.9839 0.8080 0.8677 0.9839
SEG + AV 0.8303 0.8697 0.9837 0.8267 0.8673 0.9837 0.8250 0.8653 0.9837

Table 3: Topology estimation results on INSPIRE dataset(Same parameters)

Precision Recall F


0.9416 0.9865 0.9920 0.9380 0.9863 0.9920 0.9378 0.9864 0.9920
SEG + AV 0.9450 0.9475 0.9891 0.9437 0.9469 0.9891 0.9435 0.9467 0.9891

0.8158 0.8745 0.9841 0.8146 0.8739 0.9841 0.8132 0.8723 0.9841
SEG + AV 0.8427 0.8799 0.9843 0.8391 0.8778 0.9843 0.8379 0.8762 0.9843

Table 4: Topology estimation results on INSPIRE dataset(Optimized parameters: , , , , , )

5 Discussion

Different diseases are have different pathological markers uniquely identified in a fundus image. Most of the deep learning techniques tend to only diagnose absence or presence of a single disease. This leads to extra time, effort, and extra money to make another system to diagnose another disease. In this work, we attempt to move forward to the first step to extract all the vascular features of fundus images (fig. 1). Features like Artery/Vein ratio, tortuously, bifurcation, disconnected vessels can can be extracted following this work as we have shown that we have extracted both artery/Vein mask and their topology. This work is a pipeline towards it because in our ablation studies table 1, 2, 4 we have shown that the performance of it depends on quality of priors we feed to the system. Furthermore, we have shown that with decent priors, the algorithm trained on DRIVE dataset works well with another dataset with no to minimum parameter calibration. We are confident that this work enables a solid/more principled path for automated fundus image diagnosis. In addition, in fig. 4, we can see how our flow-assignment algorithm detects possible points where vessel are broken. The size of the nodes depicts if a natural blood path has been broken or obstructed. For example, in fig. 4, we can see disconnected vessel end point get more visits than the other ones(near ONH center left of first image, top right of second and third image). We have only used the directed graph obtained from the flow-assignment algorithm, however intuitive depiction blood flow obstruction could be used in future research as a important factor in itself.

6 Conclusion and Future Work

Automated diagnosis of diseases is not fully trusted in medical field. One reason being, it not explainable enough. Thus we aimed to enable diagnose fundus image by extracting most of the features like a human.

Symbol Description
Fundus image
Vessel segmentation ground truth
Vessel segmentation likelihood (probability map)
Artery-vein ground truth
Artery-vein likelihood (prior)
Min diameter for connected components to be part of a vessel
Weights by vessel prob. map in Dijkstra’s algorithm
Weights by pixel color in Lab space in Dijkstra’s algorithm
Weights by pixel distance from background in Dijkstra’s algorithm
How much to expand a vessel path
Min number of nodes for a segment
Number of iterations of Dijkstra’s algorithm
-expansion in Dijkstra’s algorithm
False segment limit
Cost jump for prop
Optic nerve head (ONH) radius
Default branch-forward pair theta coeff.
Labeled branch-forward segments pair A/V cross-entropy coeff.
Branch cost pair scale factor
Branch cost prob. scale factor
Branchness cost scale factor
Labeled branch-forward graph cost scale factor
Destination. End nodes(nodes with degree 1) of graph
Destination. End nodes of graph with smallest threshold
Branch nodes(nodes with degree 2)
Segment(List of nodes between two branches) starting from a branch towards its neighbor
Set of segments in a graph
Adjacent Nodes of in
In neighbors of in , where is directed
Out neighbors of in , where is directed
List of (integer) neighbors in a lattice of pixel
Optic nerve center pixel
Checkpoint nodes ( + [])
List of nearest neighbors in from node
List of connected components of
Empty undirected graph
Empty directed graph
Dijkstra’s shortest path from to in a weighted graph
Weight assigned to edge (i, j) in graph
Degree of node n
Average expectation of a vessel (a path of nodes)
Artery vein propagation score
Table 5: List of symbols