1 Introduction
Retinal fundus images allow ophthalmologists to diagnose a variety of ocular and cardiovascular diseases, including diabetic retinopathy (DR) [retinal_imaging_and_analysis5660089], glaucoma [economic_analysis_pmid10868871], agerelated macular degeneration (AMD) [AMD_severity_scaleFerris2005], and the likelihood of stroke [retinal_signs_and_stroke_baker_retinal_2008, Kipli2018]. These diagnoses are based on telltale features in the fundus image that are correlated with a higher likelihood of certain diseases. For example, a higher cuptodisk ratio in the optic nerve is correlated with a higher likelihood of glaucoma [cuptodiscglaucomaBOCK2010471]. Similarly, the arteryvein ratio can be used to predict a patient’s risk for diabetic retinopathy and hypertension [retinal_imaging_and_analysis5660089, 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.031390, doi:10.1117/12.708469, Lue31, retinal_imaging_and_analysis5660089, exutates_DRJOSHI20181454, turtuosity_use_HART1999239].
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. UNet, a multilevel encoderdecoder 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 UNet as the backbone for our intermediate feature extraction tasks as in figure
1, where we use two UNets for vessel probabilitymap assignment and Artery/vein probabilitymap assignment. Usually, we use a threshold to binarize such probabilitymaps to obtain a vessel mask. Using a fixed mask can lead to a very annoying tradeoff 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 beforementioned 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 multilevelskeletonization (3.1) where we use a set of thresholds ranging from lowto 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 uniongraph. This union graph can be seen in figure 3.a,b where the cyan colored nodes are branchnodes(), golden nodes are called endnodes(). Naturally the branch nodes should only appear in the branching locations, however, we can see the uniongraph being dense. Thus, we use a novel pruning technique called uniongraphcontraction (2) to get rid of significant amount of nodes and edges keeping the overall vascular structure intact( 3.c). Then we use 2step 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 biomarker 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 UNet, directed graph topology obtained by flowassignment 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 whitecoloured 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.
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 treelike graph. Then quantify structural differences between the graph of healthy and nonhealthy patients. Fractal analysis used to characterize the extracted graphs. The paper claimed there is a significant differences between healthy and momhealthy 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 graphtheoretic 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 threestage algorithm for blood vessel segmentation. In the first stage a binary image extracted from the highpass 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 Unet 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 computeraided method to predict A/V classification. In order to reduce the differences in feature space, the Intraimage regularization and intersubject normalization applied. This method tested on DRIVE Dataset[XU20173].An attention guided Unet with atrous convolution(AAUNet) applied by [8999616] to segment the vessel and nonvessel pixels. They use AAUNet 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 vesselconstraint network(VCNet) to enhance A/V classification. This method utilizes the information of vessel distribution and edge. VCNet uses a vessleconstraint (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 backgroundprone 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 VCNet 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) arteryvein 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 )AVratio, 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 endtoend 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.
3.1.1 Extraction of LikelihoodMap
We use Convolutional Neural Network on the green channel to get vessel likelihood map as . Specifically we use a UNet architecture [Ronneberger2015UNetCN] with stochastically weighted loss [deeddyn_10.3389/fcomp.2020.00035] to train kfold 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 AVPrior
We extract likelihood for Artery/Vein (noted ) for each pixel by kfold cross validation with UNet. Such pixelwise classification is inconsistent to actually use for measures like AV ratio, bifurcation because of too many misclassification 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 multilevelskeleton(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).
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.
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:(1) 
, 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 unvisited vessels, leading to disconnected paths. In order to fix such issues we, first detect such lumps of unvisited regions, normalize edge weight based on distance to the current extracted topology(alg. 5) and rerun the algorithm in such regions only (4).
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 NPhard, thus our flow assignment algorithm 7 is robust to assign approximate blood flow.
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 branchforward method that determines what metrics to use to choose the path in a branch. One can choose to use any branchforward measure, however we choose the following:
(2) 
where is the straightness of segment when forwarding to another segment of a branch s.t. (fig. 5).
We use the branchforward 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 misalignment 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.
r
3.2.2 Tree topology cost function.
We use an intuitive cost function to optimize the topology:
(3) 
, where is the undirected 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 AA/VV in check.
(4) 
, 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 .
(5) 
, 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 outdegree should always be more than indegree. Also, when there are multiple vessels crossing, the segment that share a minimum branch forward cost () as in fig. 5 should be preferred.
The cost is incurred when the branchforward 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 branchforward method that has the A/V prior factor added along with the orientation(straightness) as mentioned in fig. 5 as follows:
(6) 
(7) 
, where s.t. . For all inneighbors 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 NPhard problem [wide_dataset:6987362], however a good directed graph from section 3.1, and arteryvein prior from UNet from section 3.1.1 makes the optimization problem much easierwe 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 3degree source and sink nodes instead of actual 4degree 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 undirected 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 subgraph 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.
(8) 
, where , and , where is just a scaling factor for adjustment if necessary.
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: ShiftUp, ShiftDown, ReverseShiftUp, and ReverseShiftDown. 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 ShiftSourceUp as shown in figure 8, and there could be another one similar but ShiftSourceUp 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 beforementioned 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 shiftup 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:ShDWNSRC (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 avpmap 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 falsepropagation 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 stoppoints: 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 4neighbor branch(crossing) end up being a 3neighbor 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 stoppoints by just adding these stopnodes 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_datasettang2011robust]. 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 AVLabel() in which each pixel is assigned 3 probabilities as a distribution—artery, vein, and background. in subscript specify likelihood prior obtained by a trained UNet, 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 groundtruth never come together in one experiment.
4.1 Transfer learning for INSPIRE dataset
For INSPIRE dataset[inspire_datasettang2011robust], 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 UNet 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 AVsegmentation groundtruth . We then train UNet 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 Nodenode comparison
is the score of nodes of graph from section 3.1 labeled based on corresponding prior in combination and mapped with groundtruth. 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 AVlabel comparison
In this approach, instead of only considering common pixels, like section 4.2.2, we compare complete pixels in the imagewe will have 3 possible labels for each pixel—artery, vein, and background. Then we use the manual AV groundtruth to calculate score. We keep record of this score before and after applying propagation algorithm.

Precision  Recall  F  

Combination 
P  P  P  R  R  R  F  F  F 
SEG + AV 
0.9181  0.9898  0.9708  0.8872  0.9895  0.9708  0.8941  0.9896  0.9708 
SEG + AV 
0.9325  0.9638  0.9684  0.9274  0.9643  0.9684  0.9287  0.9638  0.9684 
SEG + AV 
0.7814  0.8738  0.9649  0.7822  0.8723  0.9649  0.7784  0.8710  0.9649 
SEG + AV 
0.8171  0.8806  0.9653  0.8180  0.8798  0.9653  0.8143  0.8781  0.9653 
SEG + AV 
0.7780  0.8815  0.9625  0.7808  0.8795  0.9625  0.7754  0.8783  0.9625 
SEG + AV 
0.8140  0.8888  0.9631  0.8174  0.8877  0.9631  0.8126  0.8863  0.9631 


Precision  Recall  F  

Combination 
P  P  P  R  R  R  F  F  F 
SEG + AV 
0.9117  0.9947  0.9753  0.8816  0.9947  0.9753  0.8858  0.9947  0.9753 
SEG + AV 
0.9319  0.9670  0.9726  0.9299  0.9663  0.9726  0.9302  0.9665  0.9726 
SEG + AV 
0.8122  0.8850  0.9721  0.8114  0.8825  0.9721  0.8102  0.8822  0.9721 
SEG + AV 
0.8739  0.9127  0.9736  0.8724  0.9111  0.9736  0.8715  0.9107  09736. 
SEG + AV 
0.7982  0.8955  0.9681  0.7989  0.8931  0.9681  0.7967  0.8927  0.9681 
SEG + AV 
0.8305  0.9042  0.9685  0.8311  0.9031  0.9681  0.8285  0.9026  0.9685 


Precision  Recall  F  

Combination 
P  P  P  R  R  R  F  F  F 
SEG + AV 
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 
SEG + AV 
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 


Precision  Recall  F  

Combination 
P  P  P  R  R  R  F  F  F 
SEG + AV 
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 
SEG + AV 
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 

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 flowassignment 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 flowassignment 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 

Thresholds  
Fundus image  
Vessel segmentation ground truth  
Vessel segmentation likelihood (probability map)  
Arteryvein ground truth  
Arteryvein 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 branchforward pair theta coeff.  
Labeled branchforward segments pair A/V crossentropy coeff.  
Branch cost pair scale factor  
Branch cost prob. scale factor  
Branchness cost scale factor  
Labeled branchforward 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 