Advanced Image Analysis for the Assessment of Retinal Vascular Changes

FACULDADE DE ENGENHARIA DA UNIVERSIDADE DO PORTO Advanced Image Analysis for the Assessment of Retinal Vascular Changes Behdad Dashtbozorg Doctoral Program in Electrical and Computer Engineering Supervisor: Prof. Ana Maria Mendonça Co-supervisor: Prof. Aurélio Campilho May 2015 © Behdad Dashtbozorg, May 2015 This work was supported by FEDER funds through the Programa Operacional Factores de Competitividade-COMPETE and by Portuguese funds through FCT-Fundação para a Ciência e a Tecnologia in the framework of the projects PEst-C/SAU/LA0002/2011, PEst-C/SAU/LA0002/2013 and the research grant SFRH/BD/73376/2010. Advanced Image Analysis for the Assessment of Retinal Vascular Changes Behdad Dashtbozorg Doctoral Program in Electrical and Computer Engineering To my loving wife Bahareh i ii Abstract In the last decade, one of the major advances in retinal vascular imaging research has been the clear demonstration that physiological and pathological alterations in the retinal vascular network are associated with a variety of worldwide major diseases such as diabetes, hypertension and atherosclerosis. However, the clinical assessment of the retinal vascular condition is most of the times tiresome, and prone to errors, particularly if occurs in a screening environment. Recent advances in image analysis can avoid this workload and provide the ophthalmologist with objective and reproducible results useful in daily clinical practice. The retinal image analysis systems have therefore become a prominent and powerful diagnostic tools in the field of ophthalmology by detecting the changes in retinal images. It has been widely demonstrated that in diabetic retinopathy, the blood vessels often show abnormalities at early stages, as well as vessel diameter alterations. Changes in retinal blood vessels, such as significant dilatation and elongation of main arteries, veins, and their branches, are also frequently associated with hypertension and other cardiovascular pathologies. Among several characteristic signs associated with vascular changes, the Central Retinal Arteriolar Equivalent (CRAE), the Central Retinal Venular Equivalent (CRVE) and the Arteriolar-to-Venular Ratio (AVR) have been frequently used as indicators for the early detection, diagnosis, staging and follow-up of diabetes and hypertension, since they can reflect the narrowing or dilation of the retinal blood vessels. The main goal of this work is the development of an automatic system for the measurement of CRAE, CRVE, AVR and several bifurcation geometrical features. Among other image processing operations, the estimation of these features requires vessel segmentation, vessel caliber measurement, artery/vein (A/V) classification and optic disc (OD) segmentation. iii iv The CRAE, CRVE and AVR values are calculated from the calibers of the vessels inside a specific region of interest (ROI), defined as the standard ring area around the OD. As a consequence, both the localization of the optic disc center (ODC) and its diameter are required for automating the AVR calculation. For this reason, a fully automatic method based on sliding band filter is proposed; this method is able to produce useful results even in the presence of severe pathological conditions and showing a great independence from image acquisition settings. In this work, a graph-based method is proposed for the classification of retinal vessels as arteries or veins using a combination of structural information taken from the vasculature graph with intensity features from the original color image. Supervised and unsupervised techniques are introduced for the final assignment of A/V classes; the supervised approach uses linear discriminant analysis and a set of intensity features, while the unsupervised approach assigns the A/V classes using a k-means clustering algorithm and red intensity of vessel pixels. Finally, the developed approaches are integrated in RetinaCAD, Retinal Computer Aided Diagnosis system which includes vessel segmentation, vessel width estimation, optic disc segmentation and A/V classification. This system is designed to facilitate the application of retinal image analysis tools and the automated estimation of some indexes, in particular CRAE, CRVE, AVR and bifurcation geometrical features with high potential for clinical applications. Resumo Na última década, um dos maiores avanços na investigação em imagens vasculares da retina foi a clara demonstração de que alterações fisiológicas e patológicas na rede vascular estão relacionadas com uma variedade de doenças prevalentes em todo o mundo tais como diabetes, hipertensão e aterosclerose. Contudo, a avaliação clínica da condição vascular retiniana é na maioria das vezes uma tarefa cansativa e suscetível à ocorrência de erros, particularmente se a avaliação for feita num ambiente de rastreio. Avanços recentes em Análise de Imagem permitem evitar esta sobrecarga de trabalho e fornecerao oftalmologista resultados objetivos e reprodutíveis que são úteis na prática clínica diária. Os sistemas de análise de imagens da retina tornaram-se assim ferramentas de diagnóstico importantes para a deteção de alterações em imagens da retina. Na retinopatia diabética tem sido amplamente demonstrado que os vasos sanguíneos podem apresentar alterações em fases precoces, nomeadamente alterações no seu diâmetro. Alterações nos vasos sanguíneos da retina, tais como a dilatação significativa e o alongamento das principais artérias, veias e seus ramos, são também frequentemente associadas à hipertensão e outras patologias cardiovasculares. Entre vários sinais característicos associados a alterações vasculares, o Diâmatro Arterial Equivalente (CRAE), o Diâmetro Venoso Equivalente (CRVE) e o Índice Arterio-Venoso (AVR) têm sido frequentemente utilizados como indicadores para a detecção precoce, diagnóstico, avaliação do estado e seguimento de diabetes e hipertensão, uma vez que os seus valores podem refletir o estreitamento ou dilatação dos vasos sanguíneos da retina. O objetivo principal deste trabalho é o desenvolvimento de um sistema automático para a medição do CRAE, CRVE, AVR e várias características geométricas das bifurcações vasculares. Entre outras operações de processamento de imagem, a estimação destas características requer a segmentação e medição do calibre dos vasos sanguíneos, classificação dos vasos em artérias ou veias e segmentação do disco óptico. v xii CONTENTS 2.3.1 A/V Classification Evaluation . . . . . . . . . . . . . . . . . . . . . . . 24 2.4 Vascular Changes Assessment . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 2.5 ConcludingRemarks ................................ 29 3 Optic Disc Segmentation 31 3.1 VesselSegmentation ................................ 32 3.2 Optic Disc Segmentation Using Sliding Band Filter . . . . . . . . . . . . . . . . 34 3.2.1 Preprocessing................................ 35 3.2.2 SlidingBandFilter............................. 36 3.2.3 Low-resolution ODC Estimation . . . . . . . . . . . . . . . . . . . . . . 37 3.2.4 ODSegmentation.............................. 41 3.3 Results........................................ 44 3.4 ConcludingRemarks ................................ 54 4 Artery/Vein Classification 55 4.1 Graph-based A/V Classification Method . . . . . . . . . . . . . . . . . . . . . . 56 4.1.1 GraphGeneration.............................. 56 4.1.2 GraphModification............................. 59 4.1.3 GraphAnalysis............................... 62 4.1.4 A/V Class Assignment . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.2 Results........................................ 76 4.3 ConcludingRemarks ................................ 84 5 Assessment of Retinal Vascular Changes 85 5.1 Background..................................... 86 5.2 Arteriolar-to-Venular Ratio (AVR) Calculation Method . . . . . . . . . . . . . . 87 5.3 Results........................................ 95 5.4 ConcludingRemarks ................................ 102 6 Experimental Results 103 6.1 RetinaCADSystem................................. 104 6.1.1 Tools .................................... 104 CONTENTS xiii 6.1.2 Measurements ............................... 106 6.1.3 AdditionalFeatures............................. 107 6.2 EvaluationResults ................................. 108 6.2.1 Material................................... 108 6.2.2 Experimental Validation . . . . . . . . . . . . . . . . . . . . . . . . . . 108 6.2.3 ClinicalValidation ............................. 116 6.3 ConcludingRemarks ................................ 119 7 Conclusions and Future Work 121 7.1 SummaryandConclusions ............................. 121 7.2 FutureDirections .................................. 123 A Materials 125 A.1 STAREDataset................................... 125 A.2 DRIVEDataset................................... 126 A.3 VICAVRDataset .................................. 126 A.4 INSPIRE-AVRDataset............................... 126 A.5 MESSIDORDataset ................................ 127 A.6 ONHSDDataset................................... 127 B Publications 129 References 131 xiv CONTENTS List of Figures 2.1 The retina (a) Eye anatomy; (b) Retinal image. . . . . . . . . . . . . . . . . . . 8 3.1 Examples of vessel segmentation results; left column: Original images; right column: segmentation result; (a), (b) DRIVE dataset; (c), (d) STARE dataset; (e), (f) ONHSD dataset; (g), (h) INSPIRE-AVR dataset; (i), (j) MESSIDOR dataset. . . . 33 3.2 Block diagram of SBF-based optic disc segmentation. . . . . . . . . . . . . . . . 34 3.3 (a) Original image; (b) Vessel segmentation image; (c) Result of vessel pixels elimination; (f) IRG image with the initial ODC (black cross). . . . . . . . . . . 36 3.4 Schematics of the sliding band filter with 8 support region lines (dashed lines), where a simplified support region is depicted with the segmented lines. The gray region specifies a denser support region using a higher number of radial lines. . . 38 3.5 Geometric representation of camera field of view. . . . . . . . . . . . . . . . . . 39 3.6 (a) IRG image with the ROI for the low-resolution SBF (square); (b) Low-resolution SBF response on ROI; (c) 10 highest values of filter response (dots) and mean value coordinates (star); (d) Initial ODC (black cross) and new ODC candidate (graycross). .................................... 40 3.7 Geometric representation for calculating the number of support region lines. . . 41 3.8 (a) Maximum value of high-resolution SBF response (star) and the obtained band support points (dots); (b) Band support points in cropped image round OD; (c) Polar plot of smoothing result (solid line) on the band support points (dots); (d) OD boundary after smoothing on the original RGB image. . . . . . . . . . . . . 43 3.9 Comparison between proposed SBF-based method and four other methods in terms of percentage images per subjective category (ONHSD). . . . . . . . . . . . . . 47 xv xvi LIST OF FIGURES 3.10 Samples of OD segmentation in the ONHSD (solid line: results of proposed method, dots: mean of clinician boundaries). (a) Excellent (δ=0.5); (b) Good (δ=1.4); (c) Fair (δ=2.3)............................. 48 3.11 Comparison between the SBF-based method and other methods in terms of percentage of images per overlapping interval (MESSIDOR dataset). . . . . . . . . 49 3.12 Samples of OD segmentation in MESSIDOR dataset (Dashed line: results of proposed method, Solid line: manually extracted boundaries by experts). (a) S=0.97; (b) S=0.86; (c) S=0.78. ............................. 49 3.13 Samples of OD segmentation in the presence of exudates, peripapillary atrophy and blurredness; (a) DR: D3 and ME: R2; (b) DR: D3 and ME: R1; (c) DR: D3 and ME: Normal; (d) DR: D1 and ME: Normal ; (e) DR: D3 and ME: Normal; (f) DR:D3andME:R2. ................................ 51 3.14 Samples of OD segmentation in different conditions of contrast and illumination. 52 3.15 OD segmentation in MESSIDOR dataset where the initial OD detection (black cross) failed with the low-resolution ODC location (gray cross) and final ODC location (white cross); (a and b) The method failed to detect and segment the OD (S=0); (c) The method overcame the initial OD localization failure and segmented the OD correctly (S=0.90)......................... 52 3.16 Samples of OD segmentation in INSPIRE-AVR dataset (Dashed line: results of proposed method, Solid line: manually extracted boundaries by an expert). (a) S=0.97; (b) S=0.86; (c) S=0.75......................... 53 4.1 Block diagram of the proposed method for A/V classification. . . . . . . . . . . 57 4.2 Graph generation. (a) Original image; (b) Vessel segmentation image; (c) Centerline image; (d) Extracted graph. . . . . . . . . . . . . . . . . . . . . . . . . . . 58 4.3 Graph Modifications. (a), (d), (g) Typical errors; (b), (e), (e) Graph representation of the worst-case scenarios; (c), (f), (i) Final graph after modification. . . . . . . 60 4.4 (a) Original graph; (b) New graph after modification step; (c) Modified graph without vessels around the optic disc; (d) Separate subgraphs (each gray region representsasubgraph). ............................... 62 LIST OF FIGURES xvii 4.5 (a)-(c) Possible configurations for nodes of degree 2; (d)-(f) Possible configurationsfornodesofdegree3.............................. 64 4.6 (a)-(c) Possible configurations for nodes of degree 4; (d) Node of degree 5. . . . . 66 4.7 Examples of subgraph labeling where each color represent a distinct label; (a) Paired subgraph 1; (b) Paired subgraph 2; (c) Unpaired subgraph 3; (e), (e) and (f) Results of link labeling in subgraphs (a), (b) and (c). . . . . . . . . . . . . . . . 68 4.8 (a) Separate subgraphs (b) Final result of graph analysis. . . . . . . . . . . . . . 68 4.9 Examples of assigning A/V classes to the labels in each subgraph using LDA classifier; (a), (d) and (g) Graph analysis results for each subgraph; (b), (e) and (h) LDA classifier result; (c), (f) and (i) Final result of assigning A/V classes. . . . . 71 4.10 (a) LDA classifier result; (b) Final result of supervised graph-based A/V classification; (c) A/V classification result overlapped on original image (Red: arteries, Blue:veins). .................................... 71 4.11 Histograms and normal distribution fits for artery and vein pixels based on manual A/V classification in (a) Red channel; (b) Green channel; (c) Blue channel; (d) Hue; (e) Saturation; (f) Intensity. . . . . . . . . . . . . . . . . . . . . . . . . . . 72 4.12 (a) Original image; (b) Normalized red plane; (c) Intensity of vessel pixels. . . . 73 4.13 Initial cluster centroid positions on the sorted set of intensities. . . . . . . . . . . 74 4.14 Initial cluster centroid positions on the histogram of Red intensity. . . . . . . . . 74 4.15 (a) k-means clustering result (Red: artery, Blue: vein and Green: unknown); (b) Paired subgraphs; (c) Result of assigning A/V classes to paired subgraphs using k-meansalgorithm. ................................. 74 4.16 (a) Final result of unsupervised graph-based A/V classification; (b) Unsupervised A/V classification result overlapped on original image (Red: arteries, Blue: veins). 75 xviii LIST OF FIGURES 4.17 Samples of A/V classification results in INSPIRE-AVR dataset (Red: correctly classified arteries, Blue: correctly classified veins, Green: wrong classification); (a), (g) Original images; (b), (h) Manual A/V labeling; (c), (i) Supervised A/V classification results (accuracy = 97.7% and 61.2%); (d), (j) Comparison of supervised A/V classification results with manual labeling; (e), (k) Unsupervised A/V classification results; (f), (l) Comparison of unsupervised A/V classification results with manual labeling (accuracy = 89.6% and 70.1%). . . . . . . . . . . . 77 4.18 Accuracy of correctly classified vessel pixels for the supervised and unsupervised approaches in entire image and inside the ROI. . . . . . . . . . . . . . . . . . . 79 4.19 Performance of the supervised graph-based method (Red dot) and unsupervised approach (Yellow dot) compared with the results of Niemeijer’s method (Blue dot: best cut-off, Purple dot: sensitivity cut-off and Green dot: specificity cut-off); (a) INSPIRE-AVR dataset; (b) DRIVE dataset. . . . . . . . . . . . . . . . . . . 80 4.20 Samples of A/V classification results in DRIVE dataset (Red: correctly classified arteries, Blue: correctly classified veins, Green: wrong classification); (a), (g) Original images; (b), (h) Manual A/V labeling; (c), (i) Supervised A/V classification results (accuracy = 96.1% and 72.9%); (d), (j) Comparison of supervised A/V classification results with manual labeling; (e), (k) Unsupervised A/V classification results; (f), (l) Comparison of unsupervised A/V classification results with manual labeling (accuracy = 88.8% and 78.5%). . . . . . . . . . . . . . . . . . . 82 4.21 Samples of A/V classification results in CHSJ dataset (Red: correctly classified arteries, Blue: correctly classified veins, Green: wrong classification); (a), (g) Original images; (b), (h) Manual A/V labeling; (c), (i) Supervised A/V classification results (accuracy = 83.5% and 84.2%); (d), (j) Comparison of supervised A/V classification results with manual labeling; (e), (k) Unsupervised A/V classification results; (f), (l) Comparison of unsupervised A/V classification results with manual labeling (accuracy = 92.5% and 93.1%). . . . . . . . . . . . . . . . . . . 83 5.1 Block diagram of the proposed method for AVR estimation. . . . . . . . . . . . . 88 5.2 (a) Input image; (b) Binary vessel image result; (c) Distance transform result; (d) Centerlineimage................................... 89 LIST OF FIGURES xix 5.3 (a) Retinal image with the ODC (black cross) detected using the method based on the entropy of vascular directions; (b) Circular OD boundary using a fixed OD radius of 180 pixels centered on the initial ODC in (a); (c) ROI for AVR calculation (delimited by the two green circles) with a fixed OD radius; (d) OD boundary using SBF-based method; (e) Approximation of the OD boundary by a circle (radius of 215 pixels); (f) ROI for AVR calculation (delimited by the two green circles) and the estimated optic disc margin (white circle). . . . . . . . . . 91 5.4 (a) A/V classification result using the supervised graph-based method; (b) Main vessels inside the ROI (supervised AV classification and fixed OD radius); (c) A/V classification results using the unsupervised graph-based method; (d) Main vessels inside the ROI (unsupervised AV classification and OD segmentation). . . . . . 93 5.5 Region of interest divided in six concentric regions. . . . . . . . . . . . . . . . . 94 5.6 Bland-Altman plots of the agreement (a) between Observer 2 and reference; (b) between Niemeijer’s method and reference; (c) between Method 1 and reference; (d) between Method 2 and reference; (e) between Method 3 and reference. . . . . 96 5.7 Boxplot of AVR values for different methods. . . . . . . . . . . . . . . . . . . . 99 5.8 Scatter plots and regression lines between reference AVR values (a) between Observer 2 and reference; (b) between Niemeijer’s method and reference; (c) between Method 1 and reference; (d) between Method 2 and reference; (e) between between Method 3 and reference. . . . . . . . . . . . . . . . . . . . . . . . . . . 100 5.9 (a) Number of subjects with matched classification between methods and reference; (b) Number of subjects with mismatched classification between methods andreference .................................... 101 6.1 RetinaCAD graphical user interface; (a) Main screen; (b) Second display; (c) Report.105 6.2 Example images from CHSJ dataset where the RetinaCAD system fails. . . . . . 109 xx LIST OF FIGURES 6.3 Examples of obtained AVR values and AV classification results inside the ROI for different images of 2 subjects from CHSJ dataset; each column shows the results for each subject; (a), (b) results for right eye with 30° FOV; (c), (d) results for right eye with 45° FOV; (e), (f) results for left eye with 30° FOV; (g), (h) results for left eye with 45° FOV; (Green circles: ROI delimitation, White circle: estimated OD, Red: arteries, Blue: veins). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 6.4 Boxplot of AVR values for different eyes (right and left) and different FOV (45° and30°)........................................ 112 6.5 Bland-Altman plots of the agreement between AVR values from the same eye with 45° and 30° FOV for (a) right eyes; (b) left eyes. . . . . . . . . . . . . . . . . . 113 6.6 Bland-Altman plots of the agreement between AVR values from right and left eyes for (a) images with 45° FOV; (b) images with 30° FOV. . . . . . . . . . . . . . . 114 6.7 Scatter plots and regression lines for AVR values of the same eye with 45° and 30° FOV for (a) right eyes; (b) left eyes; and for AVR values from right and left eyes for (c) images with 45° FOV; (d) images with 30° FOV. . . . . . . . . . . . . . . 115 6.8 The average of AVR values with 95% confidence intervals for the subjects with different pathological status (a) Right eye; (b) Left eye; (c) Average of right and lefteyesvalues.................................... 118 List of Tables 2.1 The performance of vessel segmentation methods. . . . . . . . . . . . . . . . . . 16 2.2 The performance of optic disc segmentation methods. . . . . . . . . . . . . . . . 21 2.3 The performance of A/V classification methods. . . . . . . . . . . . . . . . . . . 25 3.1 Parameters setting defined using ONHSD. . . . . . . . . . . . . . . . . . . . . . 44 3.2 Parameters setting for low-resolution SBF and high-resolution SBF. . . . . . . . 45 3.3 Scale factors for different datasets. . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.4 Parameter settings for high-resolution SBF for MESSIDOR and INSPIRE-AVR datasets........................................ 45 3.5 Comparison between proposed SBF-based method and four other methods in terms of percentage images per subjective category (ONHSD). . . . . . . . . . . . . . 46 3.6 Comparison of the average and standard deviation of different measures between proposed method and MBM on ONHSD dataset. . . . . . . . . . . . . . . . . . 47 3.7 Comparison between the SBF-based method and other methods in terms of percentage of images per overlapping interval and average overlapping of the whole set(MESSIDORdataset). ............................. 49 3.8 Comparison of the average and standard deviation (SD) of different measures between proposed method and MBM on MESSIDOR dataset. . . . . . . . . . . . . 50 3.9 Comparison between SBF-based method and F-HLSM method in terms of percentage images per subjective category based on the ratio between MAD and estimated OD radius (MESSIDOR dataset). . . . . . . . . . . . . . . . . . . . . . . 50 3.10 The average (standard deviation) of overlapping score (S) for the images of MESSIDOR dataset with different DR and ME grades. . . . . . . . . . . . . . . . . 53 xxi 4Introduction •A fully automatic method which classifies the vessels as arteries or veins based on a graph extracted from a vascular tree. This method is able to classify the whole vascular tree and does not restrict the classification to specific regions of interest. While most of the recent methods mainly use intensity features for discriminating between arteries and veins, the proposed method uses additional information extracted from a graph which represents the vascular network. The promising results of the proposed graph-based A/V classification method on the images of three different datasets demonstrate the independence of this method with various properties, such as differences in size, quality, and camera field of view [10]. •Automatic approaches for calculating the Central Retinal Artery Equivalent (CRAE), Central Retinal Venular Equivalent (CRVE) and Arteriolar-to-Venular Ratio (AVR) in retinal images which are supported by a new approach for optic disc segmentation and by new supervised and unsupervised techniques for A/V classification. These methods are evaluated on a public dataset, where the mean error of the measured AVR values with respect to the reference was identical to the one achieved by a medical expert using a semi-automated system, thus demonstrating the reliability of the proposed methods for AVR estimation [11,12]. •RetinaCAD, a user-friendly retinal computer-aided diagnosis system, that is able to automatically detect, measure and classify two main retinal landmarks, the optic disc and the vessels. RetinaCAD can measure several vascular features that are recognized as indicators for some prevalent systemic diseases, namely CRAE, CRVE and AVR, as well as various geometrical features associated with vessel bifurcations. This application was assessed in the images of a new dataset from a local hospital, where it showed an association between AVR values and clinical information. The lower AVR values in the subjects with pathological conditions in contrast to the non-pathological ones demonstrates the potential of the system as a CAD tool for early detection and follow-up of diabetes, hypertension or cardiovascular pathologies [13]. •Publications in peer-reviewed journals and presentations in the different conferences. The list of related publications are included in Appendix B. RetinaCAD has been awarded in RecPad 2014, 20th edition of the Portuguese Conference on Pattern Recognition with the Best Poster Award and also has been selected as one of 10 finalists in iUP25k - University of Porto Business Ideas Competition. 1.3 Organizational Overview 5 1.3 Organizational Overview Together with this chapter, the research reported in this Ph.D. thesis is organized in the following chapters. Chapter 2provides a general literature review on retinal image analysis methods useful for the assessment of vascular changes, namely vessel segmentation, optic disc detection, A/V classification and assessment of vascular changes. Chapter 3introduces an automatic approach for optic disc segmentation using a multiresolution sliding band filter (SBF). Chapter 4presents an automatic artery/vein classification method based on the analysis of a graph extracted from the retinal vasculature. Final classification of a vessel segment as A/V is performed by means of supervised and unsupervised approaches. Chapter 5is devoted to the presentation of approaches for the automated assessment of vascular changes, particularly the measurement of CRAE, CRVE and AVR indexes. Chapter 6presents RetinaCAD (Retinal Computer-aided Diagnosis System) and shows the results of the proposed system evaluation as well as the clinical validation. Finally, Chapter 7concludes this thesis by providing a summary of the presented research contributions and future directions. 6Introduction Chapter 2 Retinal Image Analysis for the Assessment of Vascular Changes. Overview Retina is a light-sensitive tissue lining the inner surface of the eye (Figure 2.1(a)). When an ophthalmologist uses an ophthalmoscope to look into the eye he sees a retinal image like the one in Figure 2.1(b). In the center of the retina is the optic disc which is the beginning of the optic nerve and the entry point for the major blood vessels that supply the retina. Fovea can be seen in Figure 2.1(b) as the blood vessel-free reddish spot, in the center of the area known as the macula. The blood vessels of the retina radiate from the center of the optic disc. The walls of the retinal blood vessels are transparent and therefore the column of blood flowing in these vessels can be directly observed. The arteries appear lighter and narrower when compared to the veins. Retina is the only part of blood circulation system that can be observed directly, so any vascular changes or abnormality can be detected and can be used for screening. Some of the pathologies that affect the retina are age-related macular degeneration, glaucoma, retinopathy of prematurity, diabetic retinopathy and hypertension [14]. Some of these diseases are now amenable to automated identification and assessment. In addition, retinal blood vessel pattern can also provide information on the presence or risk of developing hypertension, diabetes, cardiovascular or cerebrovascular diseases [14]. The retinal vasculature can be observed directly and it is easily accessible to study the health of the human microcirculation in vivo. Pathological changes of the retinal vasculature, such as the appearance of microaneurysms, focal areas of arteriolar narrowing, arteriovenous nicking, and 7 8Retinal Image Analysis for the Assessment of Vascular Changes. Overview retinal hemorrhages are common fundus findings in older people, even in those without hypertension or diabetes. Recent advances in retinal image analysis have allowed reliable and precise assessment of these retinal vascular changes, as well as objective measurement of other topographic vascular characteristics such as retinal vascular widths, geometrical attributes at vessel bifurcations and vessel tracking. (a) (b) Figure 2.1: The retina (a) Eye anatomy; (b) Retinal image. In this chapter, the principles of retinal digital image analysis and the techniques for extracting the signs related to the retinal vascular changes are discussed. The methods for detecting and classifying the main retinal landmarks are critical components of circulatory blood vessel analysis systems, so it would be valuable to start with reviewing the methods for retinal vessel segmentation, optic disc segmentation and artery/vein (A/V) classification (in Sections 2.1,2.2 and 2.3). After that the methods for extracting the signs related to the retinal vascular changes are described in Section 2.4. Finally, Section 2.5 summarizes the concluding remarks. 2.1 Vessel Segmentation All retinal blood vessels originate from the optic disc and they have lower reflectance when compared to other retinal surfaces, so they appear darker than the background. Retinal vessel segmentation is one of the fundamental components of automatic retinal disease screening systems. Retinal blood vessel and their attributes, such as length, width, tortuosity and branching pattern and angles are useful for diagnosis, screening, treatment, and evaluation of various systemic diseases. 2.1 Vessel Segmentation 9 Most of the vessel segmentation methods utilise the contrast existing between the retinal blood vessel and surrounding background. However, there are several challenges in retinal vessel segmentation that makes it a non-trivial task. There are other structures in the image such as the retina boundary, the optic disc, and bright or dark spots (caused by pathologies) which can negatively influence the process of vessel segmentation. On the other hand, the narrow vessels have a very low contrast in comparison with the background which makes the discrimination even harder. In addition, the effect of central reflex in wider vessels makes it complex to distinguish one vessel from two side-by-side vessels. Vessel segmentation algorithms and methods can be divided into three main categories [15,16]: image analysis approaches, tracking-based approaches and classification-based approaches. 2.1.1 Image Analysis Approaches Image analysis approaches extract meaningful information from images and use that information for segmenting the vessels. These image analysis approaches can be organized into four categories: multi-scale approaches, centerline detection methods, region growing approaches and matched filters approaches [15,16]. Multi-scale approaches perform segmentation on different image resolutions. Thick vessels are extracted at low resolution images and then thinner structures, such as deriving branches of already segmented structures, can be segmented at higher resolution. The main advantages of this technique are the increased processing speed and the increased robustness. Qin Li et al. [17] proposed a vessel segmentation method which includes a multi-scale analytical scheme using Gabor filters and scales multiplication. Scale multiplication, which is defined as the product of Gabor filter responses at two adjacent scales, enhances the edges and filters noise. After that they use a threshold probing technique utilizing the features of the retinal vessel network which uses a line tracking algorithm to guide the selection of the threshold. Lathen et al. [18] accomplish vessel segmentation using multi-scale quadrature filtering. They combine both line and edge detection using quadrature filters across multiple scales. The filter result gives well defined vessels as linear structures, while distinct edges facilitate a robust segmentation. For quadrature filtering of 2D and 3D images, the filter kernel is applied in at least three and six uniformly distributed directions, respectively. For segmentation, the authors combine the multi-scale filtering result in a global phase map by a weighted summation, favoring scales of high strength. 10 Retinal Image Analysis for the Assessment of Vascular Changes. Overview Moghimirad et al. [19] present a multi-scale approach based on a weighted 2D medialness function. The result of the medialness function is first multiplied by the eigenvalues of the Hessian matrix in every pixel of the image in order to extract vessel’s medial-lines. After that, by extracting the centerlines of vessels and estimation of radius of vessels, the retinal vessels are segmented. Li et al. [20] used the multi-scale production of multi-scale matched filters for vessel segmentation. The scale production is used to enhance the edges and decrease noise so that some small weak vessels with low local contrast are detected with good width estimation. Saffarzadeh et al. [21] proposed a vessel segmentation technique based on a multi-scale line detection. This method uses a perceptive transform based on Weber’s law, and then reduces the impact of bright lesions by k-means clustering. Afterwards, a line operator is used in three scales for the vessels detection and the segmentation is finalized by thresholding to ignore some of the dark lesions. Centerline detection approaches extract blood vessel centerline segments then create the vessel tree by connecting these centerline segments. Different approaches are used to extract the centerline structure. Mendonça et al. [8] proposed an algorithm which starts with the extraction of vessel centerlines by using directional information from a set of four directional Differenceof-Offset-Gaussians filters. After that, a region growing process guided by some image statistics connects the candidate points. The final segmentation is obtained using an iterative region growing that integrates the contents of several binary images resulting from vessel width dependent morphological filters. Recently, Mendonça et al. extended their method for segmenting high resolution images [22]. Wu et al. [23] describe a segmentation approach for vessel centerlines based on ridge descriptors. The proposed ridge descriptor contains the normalized largest curvature and the orientations of gradients in the local neighborhood. For vessels of a certain scale, the distribution of the descriptors is assumed to have a normal distribution estimated from a training set with known ground truth. Then vessel center line segmentation can be performed based on the distance between the ridge descriptor at candidate pixels and the learned model. Region growing approaches select a set of seed points based on some predefined criteria. The initial region begins at the exact location of these seeds. The regions are then grown from these seed points to adjacent points depending on a region membership criterion like same intensity characteristics. Martinez-Perez et al. [24] presented a method based on the scale-space analysis of the first and second derivative of the intensity image which gives information about its topology and overcomes the problem of variations in contrast inherent to retinal images. The local maxima over scales of the magnitude of the gradient and the maximum principal curvature are used as 2.1 Vessel Segmentation 11 features for a region growing procedure. The growth is constrained to regions of low gradient magnitude and then the borders between regions will be defined by growing vessel and background classes without gradient restriction. Grag et al. [25] described an unsupervised, curvature-based method for segmenting the complete vessel tree from color retinal images. The vessels are modelled as trenches and the medial lines of the trenches are extracted using the curvature information derived from a curvature estimate. After that, the vessel structure is extracted using a modified region growing method, where the medial points detected by the trench detection algorithm serve as seed points. The region is grown only around a selected neighborhood of the seed point, whose size is based on the width of the largest vessel. Matched filters approaches convolve the image with multiple matched filters for extracting objects of interest, as in Sofka et al. [26] for extracting vessels in retinal images. The core of the technique is a new likelihood ratio test that combines matched filter responses, confidence and vessel boundary measures. Matched filter responses are derived in scale-space to extract vessels of widely varying widths. Vessel boundary measures and associated confidences are computed at potential vessel boundaries. The combination of these responses forms a six-dimensional measurement vector at each pixel. A training technique is used to develop a mapping of this vector to a likelihood ratio that measures the vesselness at each pixel. Finally, the new vesselness likelihood ratio is embedded into a vessel tracing framework for vessel centerline extraction. The method proposed by Ramlugun et al. [27] segments the blood vessels using an 2-D Gabor filter on a histogram-equalized image followed by hysteresis thresholding. Fathi et al. [28] presents a method based on complex continuous wavelet analysis to enhance blood vessels and to remove noise. The segmented vessels are obtained by an adaptive histogram-based thresholding procedure along with proper length filtering process. Krause et al. [29] porposed a vessel segmentation method based on the local Radon transform for the vessel smoothing. This method first enhances the contrast of blood vessels using a second-order differential operator, and then the vessels are detected by the combination of smoothing along vessel directions with contrast enhancement across them. 2.1.2 Tracking-based Approaches Tracking-based approaches apply local operators on sites known a priori to belong to a vessel and track it. These methods start from an initial point, detect vessel centerlines or boundaries 12 Retinal Image Analysis for the Assessment of Vascular Changes. Overview by analyzing the pixels orthogonal to the tracking direction. Different methods are employed for determining vessel contours or centerlines. An advantage of tracking based methods is the connectedness of vessel segments, which is not guaranteed in pixel processing based methods. Vlachos et al. [30] proposed an algorithm for vessel segmentation and network extraction in retinal images. In this method, a new multi-scale line-tracking procedure starts from a small group of pixels, derived from a brightness selection rule, and terminates when a cross-sectional profile condition becomes invalid. The multi-scale image map is derived after combining the individual image maps along scales, containing the pixels confidence to belong to a vessel. The initial vessel network is derived after map quantization of the multi-scale confidence matrix. Median filtering is applied in the initial vessel network, restoring disconnected vessel lines and eliminating noisy lines. In the final step, post-processing removes erroneous areas using directional attributes of vessels and morphological reconstruction. Niemeijer et al. [31] proposed an automated vessel linking framework that connects together separate pieces of the retinal vasculature into a connected vascular tree. To determine which vessel sections should be linked together they use a supervised cost function. Yin et al. [32] presented a probabilistic tracking method for the blood vessel detection. In the tracking process, vessel edge points are detected iteratively using local gray-level statistics and vessel continuity properties. A Gaussian-shaped curve is used for estimating the local vessel sectional intensity profiles. Local vessel structure and the edge points are obtained by means of a Bayesian method with the maximum a posteriori probability criterion. Bekkers et al. [33] proposed two vessel edge tracking algorithms which are based on an orientation score. These algorithms use invertible and non-invertible orientation scores by means of cake wavelets and Gabor wavelets. They also presented a fast method for vessel centerline tracking through a multi-scale set of noninvertible orientation scores where the multi-scale approach makes the algorithm less stable at crossings and bifurcation points. 2.1.3 Classification-based Approaches Classification-based approaches aim at finding hypotheses that explain the training data and use these hypotheses for classifying each pixel as vessel or non-vessel. Niemeijer et al. [34] extract a simple feature vector for each pixel from the green component, and then, use a k-nearest neighbor (kNN) algorithm to estimate the probability of the pixel belong to a vessel. Another supervised method, called primitive-based method, was proposed by Staal et al. [35]. This algorithm is 2.1 Vessel Segmentation 13 based on the extraction of image ridges (expected to coincide with vessel centerlines) used as primitives for describing linear segments, named line elements. Each pixel is assigned to the nearest line element to form image patches, and then classified using a set of features from the corresponding line and image patch. The feature vectors are classified using a kNN classifier. The feature selection is based on sequential forward feature selection which, in each step, selects the best feature that satisfies some criterion function and includes it in the current feature set. The set that gives the best performance is chosen. The method presented in [36] by Soares et al. also adopts supervised classification. Each image pixel is classified as vessel or non-vessel based on the pixel feature vector, which is composed of the pixel intensity and 2-D Gabor wavelet transform responses taken at multiple scales. A Gaussian-mixture model classifier (a Bayesian classifier in which each class-conditional probability density function is described as a linear combination of Gaussian functions) is then applied to obtain the final segmentation. Lupascu et al. [37] proposed a method for automated vessel segmentation in retinal images. For each pixel in the field of view of the image, a 41-D feature vector is constructed. The feature vector consists of the output of filters, vesselness, and ridgeness measures based on eigen decomposition of the Hessian computed at each image pixel, and the output of a 2-D Gabor wavelet transform taken at multiple scales. Moreover, the feature vector includes the principal curvatures, the mean curvature, and the values of principal directions of the intensity surface computed at each pixel of the green component image. The value of the root mean square gradient and the intensity within the green component at each pixel are also included in the feature vector. After that an AdaBoost classifier is trained on gold standard examples of vessel and nonvessel pixels, and then used for classifying previously unseen images. Neural networks (NN) are used for simulating biological learning and are widely used in pattern recognition mainly for classification. One of the advantages that make neural networks attractive in medical image segmentation is their ability to use nonlinear classification boundaries obtained during the training of the network. Another attractive feature of the neural nets is the ability to learn. In this regard, Marin et al. [38] presented a supervised method for blood vessel detection in digital retinal images. This method in based on a feature vector with gray-level and moment invariants-based features and uses a neural network scheme for pixel classification. Two classification stages are considered: a design stage, in which the NN configuration is decided and the NN is trained, and an application stage, in which the trained NN is used to classify each pixel as vessel or non-vessel to obtain a vessel binary image. At the end the classifier performance is 20 Retinal Image Analysis for the Assessment of Vascular Changes. Overview from hue, saturation, brightness space and also the variance in the red, green, and blue channels in a small region around the pixel. Then the most discriminant features (12 features) are selected using a sequential forward-floating search and finally each pixel is classified into rim, cup, or background using the set of selected features and a k-nearest neighbor classifier. Recently, Cheng et al. [56] proposed a method which classifies each superpixel as disc or nondisc region using histograms with enhanced contrast and texture features. Superpixels are local and coherent regions that provide local image information. The authors used the simple linear iterative clustering (SLIC) algorithm to aggregate nearby pixels into superpixel. Afterwards, superpixel classification is used for initialization of the disc boundary followed by a deformable model for getting the final contour. 2.2.5 Optic Disc Segmentation Evaluation Most of the optic disc segmentation methods are evaluated on two publicly available databases, the ONHSD and MESSIDOR datasets, which are described in Appendix A. The OD regions are manually delineated by experts for the images of these datasets. Table 2.2 shows the performance of some of mentioned OD segmentation methods. The OD segmentation performance on the MESSIDOR dataset is evaluated based on the average of overlapping score ( ¯ S) [42]. The overlapping score (S) measures the common area between the OD region obtained using the automatic method (A) and the region delimited by experts (E), being defined by S=Area(A∩E) Area(A∪E)(2.4) In ONHSD dataset four clinicians marked 24 OD boundary points and the methods are compared by the percentage of images in the excellent-fair category using a discrepancy value [46]. The discrepancy, δjon image jis defined by δj=∑ imj i−µj i σj i+ε(2.5) where µj iand σj iare, respectively, the mean and standard deviation of values obtained by four clinicians on spoke iof image j,mj iis the location of the boundary using the segmentation method on spoke iof image j, and εis a small value to prevent division by zero when the clinicians are in exact agreement and was set equal to 0.5. 2.2 Optic Disc Segmentation 21 As shown in Table 2.2, the template-based method by Aquino et al. [42] reports the highest percentage of images in the Excellent-Fair subjective category, while both Cheng et al. [56] and Giachetti et al. (2014) [45] obtained the average of overlapping score equal to 0.88 on the images of MESSIDOR dataset. Both the template-based and deformable model methods are based on the edge characteristics. The performance of these methods very much depends on the differentiation of edges from the disc and other structures. The template-based approaches are usually more robust to color, contrast and illumination variations and have better performance in the presence of severe pathological conditions when compared to the methods in other categories. Table 2.2: The performance of optic disc segmentation methods. Methods Dataset ¯ SExcellent-Fair subjective category Aquino et al. (2010) [42]ONHSD - 97% MESSIDOR - 86% Giachetti et al. (2014) [45] MESSIDOR 0.88 - "Simple" Lowel et al. (2004) [46] ONHSD - 47% "DV-Hough" Lowel et al. (2004) [46] ONHSD - 81% "Temporal Lock" Lowel et al. (2004) [46] ONHSD - 83% Yu et al. (2007) [50] MESSIDOR 0.84 - Morales et al. (2013) [54]ONHSD 0.80 - MESSIDOR 0.82 - Cheng et al. (2014) [56] MESSIDOR 0.88 - ¯ S: Average of overlapping score 22 Retinal Image Analysis for the Assessment of Vascular Changes. Overview 2.3 Artery/Vein Classification Retinal blood vessel classification into arteries and veins is a necessary phase for the automatic detection of vascular changes, and for the calculation of characteristic signs associated with systemic diseases. Several methods have been proposed for the A/V classification [57–71] based on different visual and geometrical features which are used for the discrimination between veins and arteries. Arteries are bright red while veins are darker, and in general artery calibers are smaller than vein calibers. Vessel caliber can be affected by diseases, therefore this is not a reliable feature for A/V classification. Arteries also have thicker walls, which reflect the light as a shiny central reflex strip [72]. Another characteristic of the retinal vessel tree is that, at least in the region near the optic disc (OD), veins rarely cross veins and arteries rarely cross arteries, but both types can bifurcate to narrower vessels, and veins and arteries can cross each other [72]. For this reason, tracking of arteries and veins in the vascular tree is possible, and has been used in some methods to analyze the vessel tree and classify the vessels [57–60]. A semi-automatic method for analyzing retinal vascular trees was proposed by Martinez-Perez et al. in [57], using geometrical and topological properties of single vessel segments and subtrees. First, the skeleton is extracted from the segmentation result, and significant points are detected. For the labeling, the user should point to the root segment of the tree to be tracked, and the algorithm will search for its unique terminal points and in the end, decide if the segment is an artery or a vein. Another method similar to this was proposed by Rothaus et al. [58], which describes a rulebased algorithm to propagate the vessel labels as either artery or vein throughout the vascular tree. This method uses existing vessel segmentation results, and some manually-labeled starting vessel segments. Lau et al. [59] constructed their graph over a restricted region of interest around the optic disc and then assigned the vessel labels by finding the optimal forest in the subgraphs. Joshi et al. [60] first separated their vascular graph using Dijkstra’s shortest-path algorithm to find different subgraphs. They then labeled each subgraph as either artery or vein using a fuzzy C-mean clustering algorithm. Grisan et al. [61] developed a tracking A/V classification technique that classifies the vessels only in a well-defined concentric zone around the optic disc. Then, by using the vessel structure reconstructed by tracking, the classification is propagated outside this zone, where little or no information is available to discriminate arteries from veins. This algorithm is not designed to consider the vessels in the zone all together, but rather partitions the zone into four quadrants, and 2.3 Artery/Vein Classification 23 works separately and locally on each of them. Vazquez et al. [62] described a method which combines a color-based clustering algorithm with a vessel tracking method. First the clustering approach divides the retinal image into four quadrants, then it classifies separately the vessels detected in each quadrant, and finally combines the results. Then, a tracking strategy based on a minimal path approach is applied to join the vessel segments located at different radii in order to support the classification by voting. A piecewise Gaussian model to describe the intensity distribution of vessel profiles has been proposed by Li et al. [63]. In this model, the central reflex has been considered. A minimum distance classifier based on the Mahalanobis distance was used to differentiate between the vessel types using features derived from the estimated parameters. Kondermann et al. [64] described two types of features and two classification methods, based on support vector machines and neural networks, to classify retinal vessels. One type of features is profile-based, while the other is based on the definition of a region of interest (ROI) around each centerline point. To reduce the dimensionality of the feature vectors, they used multiclass principal component analysis (PCA). Niemeijer et al. [65] proposed an automatic method for classifying retinal vessels into arteries and veins using image features and a classifier. A set of centerline features is extracted and a soft label is assigned to each centerline, indicating the likelihood of being a vein centerline pixel. Then the average of the soft labels of connected centerline pixels is assigned to each centerline pixel. They tested different classifiers and found that the k-nearest neighbor (kNN) classifier provides the best overall performance. In [66], the classification was enhanced as a step in calculating the Arteriolar-to-Venular Ratio (AVR). Zamperini et al. [67] focused on determining effective features for A/V classification. They compared color, spatial, and size features and concluded that a mix of color and position features provided the best results. By performing a greedy backward feature selection on a set of 86 features, a set of 16 features was selected. They use different linear and nonlinear classifiers where the Linear Bayes classifier has the best performance. Mirsharif et al. [68] use a threestep supervised classification method. They first enhance the retina image, then extract different pixel color features to separate major arteries from veins. Afterwards, the misclassifications at bifurcation points are corrected using information from structural characteristics of the retinal vascular tree. Muramatsu et al. [69] used linear discriminant analysis and a set of six simple features to 24 Retinal Image Analysis for the Assessment of Vascular Changes. Overview classify main vessels in a ring area around the optic disc. The feature set includes color intensities for centerline pixels and their contrast differences in a 5×5 pixel region around them inside the vessel, and in a 10×10 pixel region outside the vessel. Relan et al. [70] used a Gaussian mixture model on small vessel patches to classify the main vessels in each optic disc-centred quadrant. Recently they improved their method using a least square-support vector machine classifier [71], which classifies the main vessels in a defined region of interest around the optic disc. 2.3.1 A/V Classification Evaluation The performance of the described A/V classification methods is shown in Table 2.3 where the methods are compared based on the obtained accuracy and the area under the receiver operating characteristic (ROC) curve. Some of these methods classify the whole vascular tree while the other ones only classify the main vessels where the small vessels are removed using a threshold for the vessel calibers. It should be noted that the vessel caliber threshold for the selection of main vessels is not consistent for all methods. The mentioned A/V classification methods are evaluated on different sets of images, and on different region of interest (ROI). The ROI is defined as the ring area between the circumferences with various radii around the optic disc. The results obtained by Grisan et al. [61] were compared with those provided by a manual classification on a validation set of 443 vessels on the concentric ring between the circles with radii 2×ROD and 4×ROD, where ROD is the optic disc radius. They reached an overall classification accuracy of 87.6%, which increases to 93.3% if only the diagnostically important retinal vessels are considered. For evaluation,Vazquez et al. [62] used VICAVR dataset with 58 images and they reached a classification accuracy of 88.8% for the vessel segments in the region between circles with radii 2×ROD and 3×ROD. The overall classification accuracy in the method proposed by Li et al. [63] was 85.5%, when evaluated in 505 vessel segments. Kondermann et al. [64] showed that the neural network classifier using ROI feature vector can classify 95.32% vessel pixels correctly, but when using the output of a segmentation algorithm instead of hand-segmented images as basis for the classification the results deteriorate by 10% on average. The method proposed by Niemeijer et al. [65] was used on the 20 images of the DRIVE dataset, an area under the receiver operator characteristic curve of 0.88 for correctly assigning centerline pixels of main vessels was obtained. On the same dataset Mirsharif et al. [68] reported the accuracy values of 0.84 and 0.90 for the entire image and the ring area between to circles with 2.3 Artery/Vein Classification 25 2×ROD and 3×ROD radii, while Muramatsu et al. [69] achieved the accuracy of 0.93 for main vessels inside the same region. Table 2.3: The performance of A/V classification methods. Methodology Dataset Region Vessel type ACC AUC Grisan et al. (2003) [61] 24 images Ring area (2ROD <d<4ROD)Main 0.933 - All 0.876 - Li et al. (2003) [63] 505 segments - - 0.855 - Kondermann et al. (2007) [64] 4 images Entire retina area Main 0.953 - Niemeijer et al. (2009) [65] DRIVE Entire retina area Main - 0.88 Niemeijer et al. (2011) [66] INSPIRE-AVR Ring area (2ROD <d<3ROD) All - 0.84 Vazquez et al. (2013) [62] VICAVR Ring area (2ROD <d<3ROD) All 0.8768 0.89 Joshi et al. (2014) [60] 50 images Entire retina area Main 0.9642 - All 0.9144 Lau et al. (2013) [59] 2446 images Ring area (2ROD <d<5ROD) All 0.989 - Relan et al. (2013) [70] 35 images Ring area (2ROD <d<3ROD) Main 0.92 - Relan et al. (2014) [71] 70 images Ring area (2ROD <d<3ROD)Main 0.9488 - Ring area (2ROD <d<4.5ROD) 0.9396 Mirsharif et al. (2013) [68] DRIVE Entire retina area Main 0.8405 - Ring area (2ROD <d<3ROD) 0.9016 Muramatsu et al. (2011) [69] DRIVE Ring area (2ROD <d<3ROD) Main 0.93 - Zamperini et al. (2012) [67] 42 images Ring area (2ROD <d<3ROD) Main 0.931 - ROD: Radius of optic disc d: distance from optic disc center ACC: Accuracy AUC: Area under curve 26 Retinal Image Analysis for the Assessment of Vascular Changes. Overview 2.4 Vascular Changes Assessment Vessel dilation is a well-known phenomenon in diabetes and significant dilation and elongation of retinal arterioles, venules, and their macular branches occur in the development of diabetic macular edema that can be linked to hydrostatic pressure changes. Retinal arteriolar narrowing and vessel color changes are known responses of hypertension. These vessel changes are not limited to diabetes and hypertension, and are also associated with other cardiovascular and cerebrovascular diseases as well [2]. Many methods have been proposed for retinal vessel segmentation and classification, but relatively few have been described for fully automatic detection and analysis of changes in the vasculature. Among several sings related to the vascular changes, the Arteriolar-to-Venular Ratio (AVR), which is a parameter derived from vessel caliber measurements in a specific region of retinal images, is generally used as a descriptor of generalized arteriolar narrowing. A decreased ratio between the width of retinal arteries and veins (AVR) is well established to be predictive of stroke and other cardiovascular events in adults, as well as an indicator of diabetes and hypertension. Ruggeri et al. [73] proposed a method that starts with the detection of the vessel structure by highlighting the vessel network and then uses a tracking algorithm. After locating the optic disc (OD), the concentric zone around the optic disc is partitioned into four quadrants and in every quadrant a fuzzy C-mean classifier labels the vessel pixels as artery or vein. Then, the AVR value is obtained by estimating vessel calibers in the region around the optic disc. Tramontan et al. improved this algorithm by enhancing the vessel tracking phase and the structural artery/vein discrimination features [74]. Niemeijer et al. [66] proposed an automated method which combines tobogganing and vessel pixel classification for the measurement of vessel calibers and uses supervised position regression for the detection of the optic disc center. For classifying the vessels, the authors tested different classification approaches and found that the linear classifier has better performance. Finally AVR values are calculated using Knudtson’s revised formula [75]. The method proposed by Muramatsu et al. [69] includes optic disc and vessel segmentation, vessel classification using a linear discriminant classifier and the selection of two pairs of major vessels in the upper and lower temporal regions for AVR calculation. SIRIUS [76] (System for the Integration of Retinal Images Understanding Services) is a web application which has been developed to be used by specialists for computing the AVR value. In this semi-automatic application, the specialist selects the optic disc and three circles, concentric to 2.4 Vascular Changes Assessment 27 the optic disk. Then the retinal vessel tree is automatically segmented. Afterwards, the intersect points of the three concentric circles and the extracted retinal vessel tree are used for the estimation of vessel calibers. In the next step, the specialist manually labels each of those selected points into vein or artery, and finally the AVR value is calculated. Schuster et al. [77] developed a semi-automated image recognition and analysis tool for the determination of the AVR value in retinal images. They use thresholding for vessel recognition and vessel width is calculated in a semi-automated procedure. Bhuiyan et al. [78] presented a semiautomated evaluation tool which assists physicians in measuring vessel caliber from poor quality images. The system uses texture and edge information to measure vessel caliber. The graders have the ability to remove the parts of measured widths which are incorrect. The obtained vessel widths were used for calculating the Central Retinal Arterial Equivalent (CRAE) and Central Retinal Venular Equivalent (CRVE). More recently, a semi-automated system, the Singapore Eye Vessel Assessment (SIVA) system has been proposed for quantifying the retinal vasculature morphology [79]. This system uses automated methodologies for detecting important retinal landmarks and construct a representation of the vasculature. Then, trained graders are able to edit the detected structures, after which several measurements can be obtained. Fritzsche in [80] described a fully automated method that finds changes in retinal vessels including changes in vessel width. Many prior methods considered global properties of the retina [81,82]. These methods captured summary descriptions, such as the average width or the ratio of widths of the retinal vasculature, and compared these measures with those from the same individual obtained in previous visits or with population distributions. A new and accurate method to measure the width of retinal blood vessels in fundus photography is proposed in [83] by Xu et al. which is based on a graph-theoretic algorithm. The two boundaries of the same blood vessel are segmented simultaneously by converting the twoboundary segmentation problem into a two-slice, three-dimension surface segmentation problem, which is further converted into the problem of computing a minimum closed set in a node-weighted graph. In addition to the signs related to vessel width, there are other vascular signs which can provide useful information to aid prevention and treatment of diseases. The geometry of a vascular bifurcation can be summarized by the bifurcation angle and the bifurcation coefficient. Bifurcation angle is the angle subtended between the two daughter arterioles at the vascular junction. Changes 28 Retinal Image Analysis for the Assessment of Vascular Changes. Overview in bifurcation angle may reflect alterations in blood flow, endothelial dysfunction, and attenuation in oxygen saturation [84]. Increased angles have been associated with decreased retinal blood flow, whereas decreased angles are associated with ageing and hypertension [85]. The bifurcation coefficient measures the changes in the total cross-sectional area across the bifurcation which is the relationship between the diameter of the trunk vessel and the diameter of the two daughter vessels. An increased bifurcation coefficient represents wider branch vessels, and a decreased branching coefficient indicates narrower branch compared with the trunk vessel. Any changes in the bifurcation geometry are known to occur with increasing age and in diseased coronary arteries [86]. Normal retinal blood vessels are straight or gently curved. In some diseases, the blood vessels become tortuous, i.e. they become dilated and take on a twisting path. The dilation is caused by radial stretching of the blood vessel and the serpentine path occurs because of longitudinal stretching. The tortuosity may occur only in a small region of the retinal blood vessels, or it may involve the entire retinal vascular tree [87]. It has been shown that there is an association between arteriolar tortuosity and retinopathy of prematurity [88]. The degree of tortuosity of a vessel can be calculated as the ratio between the distance a vessel travels from two vessel points, and the shortest distance between the same points drawn by a straight line [86]. Neural network based retinal image analysis [89] and integrated analysis of vascular and nonvascular changes from color retinal fundus image sequences by Narasimha-Iyer et al. [90] are also major notable efforts in integrating the segmentation of the different retinal features. One way for analyzing the branching pattern is by using mathematical techniques, such as fractal and local fractal dimension. Fractals are based on the concept of self-similarity of spatial geometrical patterns despite a change in scale or magnification so that small parts of the pattern exhibit the pattern’s overall structure. The fractal dimension describes how thoroughly the pattern fills two-dimensional spaces [86]. Previous studies on the analysis of branching patterns demonstrated the fractal nature of the retinal blood vessel network and that the fractal dimension can be used to identify a pathology [91,92]. 2.5 Concluding Remarks 29 2.5 Concluding Remarks Assessment of retinal vascular changes can provide clinically useful information to aid prevention and management of diseases such as diabetes and hypertension. Detecting changes in individual vessels requires a number of capabilities. Algorithms are needed to automatically extract the vasculature tree, accurately locate vessel boundaries, classify arteries and veins, detect the optic disc, determine region of interest, and ultimately to calculate the signs related to vascular changes. The methods that were described in this chapter still have some limitations. Most of these methods cannot achieve proper results in images with different characteristics such as size and camera field of view. The performance of these methods decreases in the images with variations of contrast and illumination, the presence of exudates and peripapillary atrophy caused by diabetic retinopathy, risk of macular edema, and the blurredness of images due to severe cataracts. Another limitation is that the mentioned systems and methods are mostly semi-automatic and require user’s interaction. In the next chapters, new retinal image analysis techniques are proposed, namely vessel caliber estimation, artery/vein classification and optic disc segmentation. A pipeline of these methods allows the computation of important vessel related indexes, namely the Central Retinal Arteriolar Equivalent (CRAE), Central Retinal Venular Equivalent (CRVE) and Arteriolar-to-Venular Ratio (AVR), as well as various geometrical features associated with vessel bifurcations, as will be described in Chapter 5. Besides the above mentioned developments in the assessment of retinal vessels, other areas of research are also active in fundus image analysis [86,93], namely: • Optic nerve head analysis; • Arteriolar color changes (copper-wiring/sliver-wiring); • Detection of retinal lesions; • Detection of irregularly shaped hemorrhages; • Detection of rare but major pathology such as neoplasms and scarring; • Detection of lesion distribution patterns, for example drusen; • Segmentation of atrophy, including geographic atrophy; • Content-based image retrieval for abnormality detection; • Change over time detection for abnormality assessment. 36 Optic Disc Segmentation (a) (b) (c) (d) Figure 3.3: (a) Original image; (b) Vessel segmentation image; (c) Result of vessel pixels elimination; (f) IRG image with the initial ODC (black cross). 3.2.2 Sliding Band Filter The SBF is a member of the convergence index (CI) filter class, whose output is a measure of the convergence degree of the gradient vectors calculated for the pixels belonging to the filter support region [96–98]. For each pixel with spatial coordinates (x,y), the convergence index (C) is defined by C(x,y) = 1 M∑ (k,l)∈R cosθi(k,l)(3.4) where θi(k,l)is the orientation angle of the gradient vector at point (k,l)with respect to the line, with direction i, that connects (k,l)to (x,y).Mis the number of points in the filter support region R. Distinct members of CI filters use different definitions for the support region, R. The support region of the coin filter (CF) is a circle with variable radius [99], while the support region of the iris 3.2 Optic Disc Segmentation Using Sliding Band Filter 37 filter (IF) can change in each direction [100]. A ring shaped region with varying radius and fixed width is the support region of the adaptive ring filter (ARF) [101]. The support region of the SBF is a fixed width band whose position in each direction changes for maximizing the convergence index value in each point [96,97]. For all these filters, the support region is usually restricted to a set of radial lines emerging from the point where the filter is being applied to, and equally distributed over a circular region centered at this point. The SBF has a more generic formulation in comparison with other CI filters, which is desirable for OD segmentation due to the fact that the shape of OD differs from an exact rounded area [102]. The SBF can also be parameterized to use a narrow band and ignore the gradient information at the center of the OD, thus reducing the interference of vessels. The SBF response at a pixel of interest P(x,y)is defined by (3.5) and (3.6) SBF(x,y) = 1 N N−1 ∑ i=0 Cmax(i)(3.5) Cmax(i) = max Rmin≤n≤Rmax 1 d n+d ∑ m=n cosθi,m!(3.6) where Nis the number of support region lines, θi,mrepresents the angle of the gradient vector at the point mpixels away from Pin direction i,dcorresponds to the width of the band, and Rmin and Rmax represent, respectively, the inner and outer sliding limits of the band. The number of support region lines (N) controls the resolution and computational cost of the SBF. The high-resolution SBF uses a higher number of support region lines, making the computation more time consuming while increasing the sensitivity of the filter to neighbouring changes; the low-resolution SBF uses a smaller number of support region lines which reduces the computation time as well as filter accuracy. Figure 3.4 presents a schematic representation of the SBF concept where the search region is defined by the inner (Rmin) and outer (Rmax) circles and the set of dpoints corresponding to maximal convergence is shown in gray; in this example only 8 support region lines are considered. 3.2.3 Low-resolution ODC Estimation As a first step, an initial ODC location is estimated and used for defining the ROI where the SBF is applied. In this phase, both an SBF with a small number of support region lines and a downsampled image are used for reducing the computation time. The filter response gives a more accurate location for the ODC to be used in the next phase. 38 Optic Disc Segmentation Figure 3.4: Schematics of the sliding band filter with 8 support region lines (dashed lines), where a simplified support region is depicted with the segmented lines. The gray region specifies a denser support region using a higher number of radial lines. 3.2.3.1 Initial Optic Disc Center Localization Regarding the computational complexity of the SBF, the OD segmentation can be speed up by focusing the SBF on a limited region of interest (ROI). For this purpose, an initial ODC is estimated using the approach based on the entropy of vascular directions by Mendonça et al. [103], using the vessel segmented image obtained in the preprocessing phase. The initial ODC, Oinit, is the center of this ROI (Figure 3.3(d)). 3.2.3.2 Image Downsampling The SBF filter is computationally expensive and its application to a high-resolution image takes a lot of time. In order to reduce the computational burden, the SBF is used twice. The first SBF is applied on a large ROI of a downsampled image in order to obtain a coarse ODC location, whose position is used for establishing a smaller ROI on the original size image for applying the second SBF. To downsample an image to a common size, an image with the resolution of 760 ×570 pixels and 45◦camera field of view (FOV) is used as reference. The images are downsampled according to the size and FOV of the reference image using the scale factor (α) that is computed by multiplying the image size scale factor (S1) and the FOV scale factor (S2). S1is defined as the quotient between the diameter of retinal image mask (d2) and the diameter of the reference image mask (d1). 3.2 Optic Disc Segmentation Using Sliding Band Filter 39 S2is defined by S1=D1 D2, where D1and D2represent, respectively, the diagonal of the reference image plane and the diagonal of the actual image plane. From Figure 3.6(d), the diagonal field of view can be calculated using Figure 3.5: Geometric representation of camera field of view. tanφi 2=Di 2fi (3.7) where φiis the field of view, Direpresents the diagonal of the image plane, and fiis the camera focal length. Using (3.7) and assuming that the focal lengths of cameras are similar (f1≈f2), the FOV scale factor is estimated by S2=D1 D2=f1tan(φ1/2) f2tan(φ2/2) f1≈f2 →S2=tan(φ1/2) tan(φ2/2)(3.8) where φ1represents the reference image FOV and φ2is the FOV of the image being processed. Finally, the scale factor for image downsampling is calculated by α=S1S2=d2tan(φ1/2) d1tan(φ2/2)(3.9) 3.2.3.3 Low-resolution SBF The image is downsampled using α, and afterwards, the SBF is applied to the ROI. The ROI is considered as an w×wsquare window centred at Oinit. Figure 3.6(b) shows the SBF result on the ROI of Figure 3.6(a). In order to find a new candidate for ODC (Ol1), the locations of K candidate points with highest filter response are listed, Q={(xi,yi),i=1,2...,K}. These points are represented by the dots in Figure 3.6(c). The outliers of this set are excluded based on two 40 Optic Disc Segmentation (a) (b) (c) (d) Figure 3.6: (a) IRG image with the ROI for the low-resolution SBF (square); (b) Low-resolution SBF response on ROI; (c) 10 highest values of filter response (dots) and mean value coordinates (star); (d) Initial ODC (black cross) and new ODC candidate (gray cross). criteria: 1) the relation of its value with the highest value of the filter response; 2) distance to the centroid of the Kpoint set. The criteria for excluding outliers are defined by (3.10) and (3.11). SBF(xi,yi)<(1−β)max 1≤i≤K(SBF(xi,yi)) (3.10) (xi−xm)2+(yi−ym)20.5>γ(3.11) where SBF(xi,yi)is the filter response at point (xi,yi)using (3.5), and xmand ymrepresent the centroid of the Kpoints. After excluding the outliers, the new centroid of the remaining points in set Qis the new candidate for ODC (Ol1), that will be the center of the new ROI to apply the high-resolution SBF in the next phase. The new ODC is shown with a gray cross in Figure 3.6(d). 3.2 Optic Disc Segmentation Using Sliding Band Filter 41 Figure 3.7: Geometric representation for calculating the number of support region lines. 3.2.4 OD Segmentation In this phase, a high-resolution SBF is applied to the original image. High-resolution SBF is a filter with a higher number of support region lines, N. The ROI for applying the SBF is a square region around the new ODC that was calculated in the previous phase. Unlike the low-resolution version of the filter, all SBF parameters in this phase are computed using the scale factor α. 3.2.4.1 High-resolution SBF For reducing the computation time, a smaller region for the ROI is used when compared with the ROI in the previous phase. For ROI determination, a default value of whpixels for the window size is then multiplied by the scale factor obtained from (3.9). The values for SBF parameters (Rmin, Rmax and d) depend on image size, so for each particular image they are obtained multiplying the low-resolution SBF ones by α. As illustrated in Figure 3.7, the number of support region lines (N) is calculated using (3.12). sin(θ/2) = h 2αRavg θ=2π/N −−−−−→N=π sin−1(h/2αRavg)(3.12) If x1, sin−1(x)≈x, equation (3.12) can be simplified to N=2παRavg h(3.13) where Ravg is an average value for OD radius in the reference image as mentioned in 3.2.3.2, and his the distance between endpoints of the support region lines. 42 Optic Disc Segmentation 3.2.4.2 Boundary Extraction After applying the SBF to the new defined ROI, the location of the point with the maximum value of filter response, (xc,yc), is selected as the final ODC. Given the detected ODC coordinates, the OD shape can be estimated using the position of the band that maximizes the convergence index response in each direction of the support region lines. To estimate the OD shape, we start by finding the positions of the sliding band points (band support points) are found corresponding to the filter maximum response. The coordinates of the band support points (X,Y)are obtained using (3.14) and (3.15). X(θi) = xc+rmax(i)×cos(θi)(3.14) Y(θi) = yc+rmax(i)×sin(θi)(3.15) where rmax(i)is the radius in direction i, which is defined by rmax(i) = argmax Rmin≤n≤Rmax 1 d n+d ∑ m=n cosθi,m!(3.16) The band support points that represent the OD boundary are shown on the original size and cropped images, respectively, in Figure 3.8(a) and Figure 3.8(b). 3.2.4.3 Boundary Smoothing In order to finalize the OD shape, a robust local regression algorithm is used for smoothing the boundary [104]. First the distances between boundary points and the final ODC are obtained (dop(i)), and afterwards a locally weighted smoothing method is applied to the set of dop(i). The smoothing process is considered local as it smooths each value by using a subset of neighbouring data points. A robust regression weight function is defined for the data points contained within the subset, which makes the process resistant to outliers. The local regression smoothing process starts by computing the regression weights for each data point in the subset. The number of data points in the subset is set equal to N/4, and the symmetric weight function is defined by w1(j) = 1−1−2j (N/4) 3!3 ,j=1,2,3,..., N 4(3.17) 3.2 Optic Disc Segmentation Using Sliding Band Filter 43 (a) (b) (c) (d) Figure 3.8: (a) Maximum value of high-resolution SBF response (star) and the obtained band support points (dots); (b) Band support points in cropped image round OD; (c) Polar plot of smoothing result (solid line) on the band support points (dots); (d) OD boundary after smoothing on the original RGB image. where w1(j)is the weight for the jth data point within the subset, and Nis the number of boundary points. Afterwards, a second-degree polynomial is employed for the weighted linear least square regression of dop(i). For preventing the distortion introduced by outliers, the distances dop(i)are smoothed again using robust weights. The computation of robust weights requires the computation of the residuals obtained in a previous smoothing step. The robust weights are given by w2(j) =     1−(ej/6M)22,ej<6M 0,ej≥6M (3.18) where ejis the residual of the jth data point produced by the previous regression smoothing algorithm, and Mis the median of absolute deviation of the residuals. 44 Optic Disc Segmentation Table 3.1: Parameters setting defined using ONHSD. Parameter Value Description d1530 px Diameter of reference image mask φ145◦FOV of reference image Rave 40 px Average of OD radius in reference image h4px Distance between endpoints of support region lines K10 Number of selected points with highest filter response β0.03 Parameter of (3.10) for excluding outliers γ8 Parameter of (3.11) for excluding outliers px: pixel The smoothing result is shown with a solid line in Figure 3.8(c). The final boundary points are specified using the smoothed distances from the ODC. Figure 3.8(d) presents the final OD boundary overlapped on the original RGB image. 3.3 Results The automatic OD segmentation method described in the previous sections was evaluated in the images of three public datasets, ONHSD, MESSIDOR, and INSPIRE-AVR which are described in Appendix A. Using these three datasets makes it possible to analyse the performance of this new approach on images with different resolution (ranging from 760 ×570 to 2392 ×2048 pixels), acquired by various fundus cameras and distinct FOVs (30◦, 45◦). Settings The image and FOV sizes in the ONHSD were considered as reference values for parameter setting. Table 3.1 shows the values of the parameters that were established using the images of the ONHSD, and afterwards applied to the two other datasets. The second column of Table 3.2 shows the values for the parameters of the low-resolution SBF, and the third column of Table 3.2 represents the formulas to obtain the parameters for the highresolution SBF. Table 3.3 presents the scale factors for the different datasets calculated using (3.9). The SBF parameters for the MESSIDOR and INSPIRE-AVR datasets, calculated using (α), are shown in Table 3.4. 3.3 Results 45 Table 3.2: Parameters setting for low-resolution SBF and high-resolution SBF. Parameter Value (low res.) Value (high res.) Description N64 64αNumber of support region lines Rmin 20 20αInner sliding band limit Rmax 60 60αOuter sliding band limit d7 7αWidth of the band w91 αwhWindow size for ROI wh−11 Default value for window size Table 3.3: Scale factors for different datasets. Dataset Image size FOV Mask diameter S1S2Scale factor (α) ONHSD 760×570 45◦530 1 1 1 MESSIDOR 1440×960 45◦900 1.7 1 1.7 2240×1488 45◦1357 2.6 1 2.6 2304×1536 45◦1440 2.7 1 2.7 INSPIRE-AVR 2392×2048 30◦2045 3.9 1.5 5.8 Table 3.4: Parameter settings for high-resolution SBF for MESSIDOR and INSPIRE-AVR datasets. Dataset Image size αN Rmin Rmax d w subset size MESSIDOR 1440×960 1.7 109 34 102 12 19 27 2240×1488 2.6 166 52 156 18 29 41 2304×1536 2.7 173 54 162 19 31 43 INSPIRE-AVR 2392×2048 5.8 371 116 348 41 63 92 ONHSD The original ONHSD dataset has 99 images. Similar to [46] and [42], the images with no discernible OD or with severe enough cataracts to prevent meaningful segmentation are excluded, leaving a final set of 90 images for assessing the SBF method. The size of the images in this dataset is equal to the reference size. Therefore the scale factor for ONHSD images is equal to 1 and there is no need to apply the high-resolution SBF. For this reason, after getting the result from the low-resolution SBF, the boundary was extracted as mentioned in 3.2.4.2. After finalizing the OD shape using the smoothing algorithm, the results were compared with those of the Circular 52 Optic Disc Segmentation (a) (b) (c) (d) Figure 3.14: Samples of OD segmentation in different conditions of contrast and illumination. (a) (b) (c) Figure 3.15: OD segmentation in MESSIDOR dataset where the initial OD detection (black cross) failed with the low-resolution ODC location (gray cross) and final ODC location (white cross); (a and b) The method failed to detect and segment the OD (S=0); (c) The method overcame the initial OD localization failure and segmented the OD correctly (S=0.90). success rate is similar to the one for the MEF method and is higher than the ones for F-HLSM and CHT methods which are 99.08% and 98.83%, respectively. The proposed algorithm was implemented in MATLAB on an Intel CPU i7-2600k, 3.40 GHz, 8 GB RAM computer. The average running time was 10.6 s per image in the MESSIDOR dataset, 3.3 Results 53 Table 3.10: The average (standard deviation) of overlapping score (S) for the images of MESSIDOR dataset with different DR and ME grades. Diabetic retinopathy Risk of macular edema All Normal D1 D2 D3 Normal R1 R2 Number 1200 540 153 247 260 971 75 154 S0.8859 0.8885 0.8865 0.8805 0.8850 0.8853 0.8776 0.8939 (0.0818) (0.0788) (0.0633) (0.0907) (0.0835) (0.0801) (0.0945) (0.0765) (a) (b) (c) Figure 3.16: Samples of OD segmentation in INSPIRE-AVR dataset (Dashed line: results of proposed method, Solid line: manually extracted boundaries by an expert). (a) S=0.97; (b) S=0.86; (c) S=0.75. Table 3.11: The average and standard deviation of different measures on INSPIRE-AVR dataset. ¯ SDC ACC TPF FPF MCC MAD/rOD Average 0.8505 0.9168 0.9958 0.9144 0.0020 0.9163 0.0897 Standard deviation 0.0864 0.0527 0.0030 0.0592 0.0025 0.0526 0.0600 while the running times for the CHT, MEF, SPC and F-HLSM methods on the same dataset were 7.36 s, 8 s, 10.9 s and 11.3 s, respectively. INSPIRE-AVR dataset In order to evaluate the performance of this new approach on high resolution images and different FOV (30◦), the proposed method was also evaluated on the 40 images of INSPIRE-AVR dataset. For segmenting the OD, the parameters of the high-resolution SBF were set based on the last row of Table 3.4. The performance of the SBF-based method is in Table 3.11. An average overlapping score of 85% was achieved for the whole dataset, while the mean MAD value was less than 9% of the OD radius. Some images depicting the manual ground-truth and the result of the SBF-based segmentation are shown in Figure 3.16. 54 Optic Disc Segmentation 3.4 Concluding Remarks The segmentation of the optic disc in retinal images is essential for the automated assessment of vascular changes. In previous sections, a new automatic methodology for OD segmentation is described which is distinct from previous approaches. This method uses the SBF in two different phases. In the first one, a low-resolution SBF is applied to the downsampled images in order to obtain an initial estimation of the ODC location, whose position is used for establishing the ROI for the high-resolution SBF in the following phase. The parameters of high-resolution SBF are adapted to the image size and camera field of view. The maximum response of the SBF gives the band support points that are used for an initial delineation of the OD boundary, which is afterwards smoothed using robust local regression. The proposed method outperforms recently published approaches for OD segmentation. The promising results on the images of three different datasets prove the independence of this approach from changes in image characteristics such as size, quality and camera field-of-view. Chapter 4 Artery/Vein Classification The classification of retinal vessels into artery/vein (A/V) is an important phase for automating the detection of vascular changes, and for the calculation of characteristic signs associated with several systemic diseases such as diabetes, hypertension, and other cardiovascular conditions. Several works on vessel classification have been proposed [57–71], but automated classification of retinal vessels into arteries and veins has received limited attention, and is still an open task in the retinal image analysis field. This chapter is mostly based in the paper "An Automatic Graph-Based Approach for Artery/Vein Classification in Retinal Images" [10]. Here, a graph-based method for automatic A/V classification is proposed. The graph extracted from the segmented retinal vasculature is analyzed to decide on the type of intersection points (graph nodes), and afterwards one of two labels is assigned to each vessel segment (graph links). Finally, intensity features of the vessel segments are measured for assigning the final artery/vein class using two different supervised and unsupervised approaches. This chapter is organized as follows. In Section 4.1 the new methods for artery/vein classification are described. The results of tests on the images of three different datasets are presented in Section 4.2, where a comparison with the manual classifications is also included. Finally, Section 4.3 summarizes the conclusions of this chapter. 55 56 Artery/Vein Classification 4.1 Graph-based A/V Classification Method Most of the methods described in Section 2.3 use intensity features to discriminate between arteries and veins. Due to the acquisition process, very often the retinal images are non-uniformly illuminated and exhibit local luminosity and contrast variability, which can affect the performance of intensity-based A/V classification methods. For this reason, we propose a method, that uses additional structural information extracted from a graph representation of the vascular network. This method follows a graph-based approach, which focuses on a characteristic of the retinal vessel tree that, at least in the region near the optic disc, veins rarely cross veins and arteries rarely cross arteries. This assumption demands for the detection of different types of intersection points, namely: bifurcation, crossing, meeting, and connecting points. A bifurcation point is an intersection point where a vessel bifurcates to narrower vessels. In a crossing point a vein and an artery cross each other. In a meeting point the two types of vessels meet each other without crossing, while a connecting point connects different parts of the same vessel. The decision on the type of the intersection points is made based on the geometrical analysis of the graph representation of the vascular structure. Figure 4.1 depicts the block diagram of the method proposed for A/V classification. The main phases are: 1) graph generation; 2) graph analysis; and 3) vessel classification. The graph generator extracts a graph from the vascular tree, and afterwards makes a decision on the type of each intersection point (graph node). Based on the node types in each separate subgraph, all vessel segments (graph links) that belong to a particular vessel are identified and then labeled using two distinct labels. Finally, the A/V classes are assigned to the subgraph labels using supervised and unsupervised techniques. The supervised approach involves a set of features and a linear classifier, while the unsupervised approach is based on the red intensity and on a k-means clustering algorithm. In the following, each one of the phases are detailed. 4.1.1 Graph Generation A graph is a representation of the vascular network, where each node denotes an intersection point in the vascular tree, and each link corresponds to a vessel segment between two intersection points. The graph generator is a three-step algorithm. First, the segmented image is used to obtain the vessel centerlines, then the graph is generated from the centerline image, and finally some additional modifications are applied to the graph. 4.1 Graph-based A/V Classification Method 57 Figure 4.1: Block diagram of the proposed method for A/V classification. 4.1.1.1 Vessel Segmentation The vessel segmentation was described in Section 3.1. The resulting segmented image is used for extracting the graph and also for estimating vessel calibers. Figure 4.2(a) illustrates the result of vessel segmentation. 4.1.1.2 Vessel Centerline Extraction The centerline image is obtained by applying an iterative thinning algorithm described in [106] to the vessel segmentation result. This algorithm removes border pixels until the object shrinks to a minimally connected stroke. Thinning of an image Iby a sequence of structuring elements {B}={B1,B2,B3,...,Bn}, can be defined by terms of the hit-and-miss transform (A⊗B=A− (A~B) = A∩(A~B)c), where Biis a 45° rotated version of Bi−1. The process starts thinning I with B1, then thinning the result with B2, and so on until Bn, and the entire process repeats until no further changes occur. The vessel centerlines from the segmented image of Figure 4.2(b) are shown in Figure 4.2(c). 58 Artery/Vein Classification (a) (b) (c) (d) Figure 4.2: Graph generation. (a) Original image; (b) Vessel segmentation image; (c) Centerline image; (d) Extracted graph. 4.1.1.3 Graph Extraction In the next step, the graph nodes are extracted from the centerline image by finding the intersection points (pixels with more than two neighbors) and the endpoints or terminal points (pixels with just one neighbor). In order to find the links between nodes (vessel segments), all the intersection points and their neighbors are removed from the centerline image and as result an image is obtained with separate components which are the vessel segments. Next, each vessel segment is represented by a link between two nodes. Figure 4.2(d) shows the graph obtained from the centerline image of Figure 4.2(c). The graph contains nodes, and at each node several links can be connected. On the other hand, any given link can only connect two nodes. Table 4.1 shows the graph notations for nodes and links which will be used in the rest of this document. The degree of a node is the number of adjacent nodes. Two nodes in a graph are called adjacent if they are connected by one link. The angle between links is defined as the magnitude of the smallest rotation that projects one of the links onto the other by considering the common node between them as the vertex. A vessel caliber 4.1 Graph-based A/V Classification Method 59 Table 4.1: Graph Notations. Notation Description NNumber of nodes in the graph ni,1≤i<NNode i DniDegree of node iwhich is the number of adjacent nodes lip,1≤p<Dnipth link of node i di j Distance between node iand j TniType of node i ∠lipliq Angle between pth and qth links of node i Wlip Vessel caliber assigned to pthlink of node i EniNumber of degree 1 nodes adjacent to node i is assigned to each link, as the average of the calibers along the corresponding vessel segment. 4.1.2 Graph Modification The extracted graph may include some misrepresentation of the vascular structure as a result of the segmentation and centerline extraction processes. As defined in [58], the typical errors are (1) the splitting of one node into two nodes; (2) a missing link on one side of a node; (3) false link. The extracted graph should be modified when one of these errors is identified. Figure 4.3 illustrates the three typical errors before and after modification. The typical errors are depicted in the left column; graphs in the middle column are the representations of the worstcase scenarios and graphs on the right side are the modified ones. Node splitting When extracting the centerline pixels in a single intersection, there are two graph nodes instead of only one. This situation, illustrated in Figure 4.3(a), creates false nodes which affect the correctness of the final result of the graph analysis phase. For addressing this problem an adaptive parameter is defined, the threshold Tns, which is used as the criterion for merging two neighborhood nodes. Tns depends on the local vessel calibers and angles according to equations (4.1) to (4.4). Tns =1 sin(α)d12+d22+2d1d2cos(α)1/2(4.1) α=min(∠li1li2,∠lj1lj2)(4.2) 60 Artery/Vein Classification ni nj Wli1 Wlj2 (a) ni nj ס݈௝ଵ݈௝ଶ lj1 lj2 li1 lj3 li2 li3 ס݈௜ଵ݈௜ଶ (b) ni li4 li3 li1 li2 (c) Node splitting ni nj Wlj3 (d) ni nj lj1 lj3 li1 lj2 (e) ni nj lj1 lj3 li1 lj2 li2 (f) Missing link ni nj dij Wlj2 (g) ni nj lj1 lj2 li1 li2 lj3 li3 (h) ni nj lj1 lj2 li1 li2 (i) False link Figure 4.3: Graph Modifications. (a), (d), (g) Typical errors; (b), (e), (e) Graph representation of the worst-case scenarios; (c), (f), (i) Final graph after modification. d1=maxWli1,Wlj2(4.3) d2=maxWli2,Wlj1(4.4) In adjacent nodes with degree 3, Figure 4.3(b), if the distance between node niand njis smaller than the threshold (di j <Tns) and if a link in one node (common link excluded) has the same orientation of another link in the other node, and the same happens with the two remaining links, then the two corresponding nodes should be merged (Figure 4.3(c)). 4.1 Graph-based A/V Classification Method 61 Missing link For solving the missing link cases (Figure 4.3(d)), the distance from a degree 1 node (endpoints) to other nodes is calculated. If this distance is less than a threshold Tml, then the nodes will be connected with a new link as shown in Figure 4.3(f).Tml estimation is based on the widths of the intervening vessels (Figure 4.3(e)), and is given by Tml =Wli1+max p∈{1,2,3}Wljp (4.5) False link Figure 4.3(g) illustrates the last situation, corresponding to an incorrect detection of a link between two nodes. This happens when two vessels are very close to each other but they do not cross, and two close nodes (di j <Tf l) are artificially created. Equation (4.6) represents the threshold Tf l for this case, which is obtained from the maximum distance of nodes in the worst-case scenario (Figure 4.3(h)). Tf l =max p∈{1,2}Wlip sin(∠lipli3)+max q∈{1,2}Wljq sin(∠ljqlj3)(4.6) For solving this case, the angle between the links connected to each node is checked. If, for at least one node, two of its links have an identical orientation and are more or less perpendicular to the third link (the common link between two nodes), then the common link is a falsely-detected link and should be removed (Figure 4.3(i)). Algorithm 1, in page 62, states the conditions for detecting errors in the graph and for performing the necessary modifications. This algorithm repeats until no further changes occur. In order to reduce the complexity of the subsequent graph analysis, all endpoints with very short links are removed. Before initiating the graph analysis phase, all vessels around the OD are removed. The optic disc area usually contains many vessels and the graph in that area is not reliable. As these vessels are not relevant for the A/V classification process, they can be removed. The result is a graph formed by several non-connected separate subgraphs. The operations in the graph modification step are illustrated in Figure 4.4. 68 Artery/Vein Classification (a) (b) (c) (d) (e) (f) Figure 4.7: Examples of subgraph labeling where each color represent a distinct label; (a) Paired subgraph 1; (b) Paired subgraph 2; (c) Unpaired subgraph 3; (e), (e) and (f) Results of link labeling in subgraphs (a), (b) and (c). (a) (b) Figure 4.8: (a) Separate subgraphs (b) Final result of graph analysis. 4.1 Graph-based A/V Classification Method 69 Table 4.3: List of features measured for each centerline pixel. Nr. Features 1-3 Red, Green and Blue intensities of the centerline pixels. 4-6 Hue, Saturation and Intensity of the centerline pixels. 7-9 Mean of Red, Green and Blue intensities in the vessel. 10-12 Mean of Hue, Saturation and Intensity in the vessel. 13-15 Standard deviation of Red, Green and Blue intensities in the vessel. 16-18 Standard deviation Hue, Saturation and Intensity in the vessel. 19-22 Maximum and minimum of Red and Green intensities in the vessel. 23-30 Intensity of the centerline pixel in a Gaussian blurred (σ=2,4,8,16) of Red and Green plane. 4.1.4.1 Supervised A/V Class Assignment Approach As a result of the acquisition process, very often the retinal images are non-uniformly illuminated and exhibit local luminosity and contrast variability. In order to make the classifier more robust, each image is processed using the method proposed by Foracchia et al. [107], which normalizes both luminosity and contrast based on a model of the observed image. Luminosity and contrast variability in the background are estimated and then used for normalizing the whole image. For each centerline pixel, the 30 features listed in Table 4.3 are measured and normalized to zero mean and unit standard deviation. Some of these features were used previously in [65, 66]. We have tested the most commonly used classifiers, namely linear discriminant analysis (LDA), quadratic discriminant analysis (QDA), and k-nearest neighbor (kNN), on the INSPIREAVR dataset. In both LDA and QDA classifiers, it is assumed that the conditional probability density functions of both classes are normally distributed and the prior probabilities are equal. For feature selection, we have used sequential forward floating selection [108], which starts with an empty feature set and adds or removes features when this increases or decreases the performance of the classifier using the classification accuracy as criterion. Table 4.4 shows the performance of intensity-based classifiers on the INSPIRE-AVR image dataset using 2-fold cross-validation. The optimal value of parameter k(number of nearest neighbors) in KNN classier is obtained equal to 5. This table contains the accuracy values obtained when these classifiers are used for centerline pixel classification, and also for labeling vessel segments (links). The LDA classifier provided the best results and was selected for the A/V classification phase, using the set of selected features (1-2, 7, 10, 12-14, 16-17, 19-20 and 23-30). The trained classifier is used for assigning the A/V classes to each one of the subgraph labels. First, each centerline pixel is classified into A or V classes (Figure 4.10(a)), then for each label 70 Artery/Vein Classification Table 4.4: Performance evaluation and comparison of individual intensity-based classifiers (INSPIRE-AVR dataset). Classifier Selected features Most important feature Accuracycenterline pixels Accuracylabeling links LDA No. 1, 2, 7, 10, 12-14, 16, 17, 19, 20, 23-30 No. 1 75.2% 79.7% QDA No. 1-3,7,8,13,14,19,21-23 No. 1 73.2% 75.9% kNN No. 1-3,7-9,13,14,16,20,23 No. 13 68.3% 73.0.% (Ci j,j=1,2) in subgraph i, the probability of its being an artery is calculated based on the number of associated centerline pixels classified by LDA to be an artery or a vein. The probability of label Ci jto be an artery is Pa(Ci j)=na Ci j /(na Ci j +nv Ci j ) where na Ci j is the number of centerline pixels of a label classified as an artery and nv Ci j is the number of centerline pixels classified as a vein. For each pair of labels in each subgraph, the label with higher artery probability will be assigned as an artery class, and the other as a vein class. Finally, to prevent a wrong classification as a result of a wrong graph analysis, we calculate the probability of being an artery or a vein for each link individually. The probability of a link (li) being an artery (Pa(li)) is computed as Pa(li) = na li/na li+nv li, and the probability of being a vein ((Pv(li)) is computed as Pv(li) = nv li/na li+nv li, where na liis the number of centerline pixels of link (li) classified as an artery and nv liis the number of centerline pixels classified as a vein. If the probability of being an artery is higher than 0.9 (Pa(li)≥0.9) then the link will be assigned as an artery, and if Pv(li)≥0.9 then it will be assigned as a vein, without considering the result of the graph analysis. Figure 4.9 shows some examples of using the LDA result for the assignment of A/V classes to each subgraph. The A/V classification result on the link inside the ellipse in Figures 4.9(a)-4.9(c) emphasizes the impact of graph analysis on the final results of A/V classification. As it can be seen in Figure 4.9(b), if the result of LDA classifier will be used only, then the link inside the ellipse will be classified wrongly as artery, but the graph-based method correctly classified this link as vein using the labels obtained in the graph analysis step. The final result of assigning a class to the link centerline pixels is shown in Figures 4.10(b) and 4.10(c). The red color represents arteries and blue color represents the veins. 4.1 Graph-based A/V Classification Method 71 (a) (b) (c) (d) (e) (f) (g) (h) (i) Figure 4.9: Examples of assigning A/V classes to the labels in each subgraph using LDA classifier; (a), (d) and (g) Graph analysis results for each subgraph; (b), (e) and (h) LDA classifier result; (c), (f) and (i) Final result of assigning A/V classes. (a) (b) (c) Figure 4.10: (a) LDA classifier result; (b) Final result of supervised graph-based A/V classification; (c) A/V classification result overlapped on original image (Red: arteries, Blue: veins). 72 Artery/Vein Classification 4.1.4.2 Unsupervised A/V Class Assignment Approach Since images of different datasets have diverse properties, the LDA classifier in the supervised A/V class assignment approach requires a computationally demanding training phase for each dataset. This requirement prevents the method from achieving the expected A/V classification performance when the classifier is trained with images from a different set. In order to overcome this limitation, a new unsupervised approach for the final A/V class assignment is developed using k-means clustering algorithm. Figure 4.11 shows the histogram and normal distribution fit for artery pixels and vein pixels in different color planes using the manual A/V classification. As illustrated in this figure, the histogram of red intensity shows the best discrimination between artery pixels and vein pixels, although with some overlap. For this reason, we selected the red component for using in the k-means clustering algorithm. (a) (b) (c) (d) (e) (f) Figure 4.11: Histograms and normal distribution fits for artery and vein pixels based on manual A/V classification in (a) Red channel; (b) Green channel; (c) Blue channel; (d) Hue; (e) Saturation; (f) Intensity. 4.1 Graph-based A/V Classification Method 73 First from the original color image (Figure 4.12(a)), a normalized intensity image from red plane is obtained (Figure 4.12(b)) and the red intensity for all vessel pixels are extracted (Figure 4.12(c)) and stored in a set, I. The elements of the obtained set are sorted in ascending order which is used for determining three cluster centroids, Cv,Cuand Ca, that allow the initializing of ak-means algorithm for clustering each vessel pixel into one of three classes cases: 1) Artery; 2) Vein; 3) Unknown. (a) (b) (c) Figure 4.12: (a) Original image; (b) Normalized red plane; (c) Intensity of vessel pixels. As retina arteries normally appear thinner and brighter red than the corresponding veins with a normal artery-to-vein caliber ratio of 2:3 [109], to compute the initial centroids the sorted set of intensities is divided into 7 intervals, each one containing the same number of pixels. The first 3 intervals are considered as initial vein cluster, the 2 last intervals belong to initial artery cluster and the 2 middle intervals are initially considered as unknown cluster. The different number of intervals in the artery and vein classes derives from the fact that veins are larger than arteries, so we have more vein pixels than artery pixels. All intensities in the 2 middle intervals are associated with the unknown class for the case of uncertainty. As it is illustrated in Figure 4.13, the initial centroids, Cv,Cuand Ca, are set equal to the center of vein, unknown and artery initial clusters, respectively. Figure 4.14 shows the obtained initial cluster centroid positions on the histogram of Red intensity. Using the k-means algorithm and the obtained initial centroids all vessel pixels are clustered as artery, vein or unknown as shown in Figure 4.15(a). Then, the probability of a label being an artery is calculated based on the relation between the number of pixels in each cluster. Subsequently, in each paired subgraph (Figure 4.15(b)), the label with higher artery probability will be considered as an artery, and the other one as a vein, where the result is illustrated in Figure 4.15(c). 74 Artery/Vein Classification Figure 4.13: Initial cluster centroid positions on the sorted set of intensities. Figure 4.14: Initial cluster centroid positions on the histogram of Red intensity. (a) (b) (c) Figure 4.15: (a) k-means clustering result (Red: artery, Blue: vein and Green: unknown); (b) Paired subgraphs; (c) Result of assigning A/V classes to paired subgraphs using k-means algorithm. 4.1 Graph-based A/V Classification Method 75 In the next step, the two thresholds are recalculated based on the result of A/V assignment in paired-subgraphs. The threshold values for arteries (Ta) and veins (Tv) are set as Ta=µa−σa(4.7) Tv=µv+σv(4.8) where µais the average intensity and σais the standard deviation of all pixels in the artery subgraphs previously classified, respectively, and µvand σvhave identical definition for the vein subgraphs. Afterwards, the classification process is repeated for all vessel pixels based on the obtained threshold. Each pixel (p) with intensity of Ipis classified as following: For each pixel (p)         if Ip≤Tv⇒p∈Veinclass if Ip≥Ta⇒p∈Arteryclass if Tv<Ip<Ta⇒p∈Unknownclass (4.9) For the subgraphs, each vessel pixel is counted as a vein or an artery using the threshold values and the probability of each label to be an artery is calculated. Then for each label in each unpaired subgraph if the probability of being artery is higher than 0.5 then the label will be assigned as artery or otherwise it will be assigned as vein; and for each pair of labels in paired subgraphs, the label with higher artery probability is assigned as an artery class, and the other as a vein class. Finally, to prevent a wrong classification of a link as a result of an error in the analysis of the graph, A/V probability for each individual vessel (graph link) is also calculated. If one of these probabilities is higher than 0.8, the vessel is considered as belonging to that class independently of the result derived from the subgraph classification procedure. Final result of unsupervised A/V classification is shown in Figure 4.16. (a) (b) Figure 4.16: (a) Final result of unsupervised graph-based A/V classification; (b) Unsupervised A/V classification result overlapped on original image (Red: arteries, Blue: veins). 76 Artery/Vein Classification 4.2 Results The automatic methods described in the previous sections were tested on the images of two datasets, DRIVE [110] and INSPIRE-AVR [111] which are described in Appendix A. A manual artery/vein labeling was performed by an expert on the 20 images of the DRIVE test set and for the 40 images of the INSPIRE dataset. In following, the proposed methods were evaluated on these datasets and finally the performance of both supervised and unsupervised A/V classification approaches were compared using 25 images of a dataset from a local hospital (CHSJ dataset). INSPIRE-AVR dataset For applying the supervised graph-based method on the images of INSPIRE-AVR dataset, 2-fold cross-validation is used. The images are randomly assigned to two sets S1and S2, so that both sets were of equal size. Then the classifier is trained on S1and it is tested on S2, followed by training on S2and testing on S1. For training the LDA classifier, 15,000 labeled centerline pixels were randomly selected from each set. Some results of the proposed A/V classification methods on this dataset are shown in Figure 4.17. Table 4.5 shows the performance evaluation of the individual LDA classifier for classifying vessel segments in INSPIRE-AVR dataset, and the results obtained using the combination of graph-based classification with LDA and unsupervised graph based approach. The analysis of these values shows that both supervised and unsupervised graph-based methods outperform the accuracy of the LDA classifier alone. For further evaluation of the graph-based method, we consider a semi-automatic approach by manually assigning the A/V classes to the labels of each subgraph. The results obtained are included in the first row of Table 4.5. For evaluating the proposed methods, the accuracy is calculated both for centerline pixel classification and for vessel pixel classification. Table 4.6 and Table 4.7 show the accuracy values for centerline and vessel pixels in the entire image, as well as for the pixels inside the region of interest (ROI) for supervised and unsupervised approaches, respectively. The ROI is usually defined for the calculation of the arteriolar-to-venular ratio; the ROI is the standard ring area within 0.5 to 1.0 disc diameters from the optic disc margin [75]. Each row of these tables contains the accuracy values calculated using different ranges for vessel calibers. The first rows contain the results for all the vessels, while the remaining rows present the results for vessels with caliber higher than 5, 10, 15 and 20 pixels. 4.2 Results 77 (a) (b) (c) (d) (e) (f) (g) (h) (i) (j) (k) (l) Figure 4.17: Samples of A/V classification results in INSPIRE-AVR dataset (Red: correctly classified arteries, Blue: correctly classified veins, Green: wrong classification); (a), (g) Original images; (b), (h) Manual A/V labeling; (c), (i) Supervised A/V classification results (accuracy = 97.7% and 61.2%); (d), (j) Comparison of supervised A/V classification results with manual labeling; (e), (k) Unsupervised A/V classification results; (f), (l) Comparison of unsupervised A/V classification results with manual labeling (accuracy = 89.6% and 70.1%). 84 Artery/Vein Classification 4.3 Concluding Remarks The classification of arteries and veins in retinal images is essential for the automated assessment of vascular changes. In the previous sections, new automatic methodologies are described for the classification of retinal vessels into arteries and veins which is distinct from previous publications. One major difference is the fact that the proposed methods are able to classify the whole vascular tree and does not restrict the classification to specific regions of interest, normally around the optic disc. While most of the previous methods mainly use intensity features for discriminating between arteries and veins, the proposed methods use additional structural information extracted from the graph representation of the vascular network. The information about node degree, the orientation of each link, the angles between links, and the vessel caliber related to each link are used for analyzing the graph, and then decisions on type of nodes are made (bifurcation, crossing, or meeting points). Next, based on the node types, the links that belong to a particular vessel are detected, and finally A/V classes are assigned to each one of these vessels using supervised and unsupervised approaches. The supervised graph-based method with LDA and the unsupervised graph-based method using k-means clustering outperform the accuracy of the LDA classifier using intensity features, which shows the relevance of using structural information for A/V classification. Furthermore, the proposed approaches in comparison with other recently punished methods, achieved better results. The promising results of proposed graph-based A/V classification methods on the images of three different datasets demonstrate the independence of these methods in A/V classification of retinal images with different properties, such as differences in size, quality, and camera angle. On the other hand, the high accuracy achieved by these methods, especially for the largest arteries and veins, confirm that this A/V classification methodology is reliable for the calculation of several characteristic signs associated with vascular alterations. Chapter 5 Assessment of Retinal Vascular Changes Automated detection of retinopathy in eye fundus images using digital image analysis has huge potential benefits, allowing the examination of a large number of images in less time, with lower cost and reduced subjectivity than current observer-based techniques. Another advantage is the possibility to perform automated screening for pathological conditions, such as diabetic retinopathy, in order to reduce the workload required of trained manual graders [86]. The Arteriolar-to-Venular Ratio (AVR) is an index used for the early diagnosis of diseases such as diabetes, hypertension or cardiovascular pathologies. Here, we present three automatic approaches for the estimation of the AVR in retinal images that result from the combination of different methodologies in some of the processing phases used for AVR estimation. Each one of these methods includes vessel segmentation, vessel caliber estimation, optic disc detection or segmentation, region of interest determination, vessel classification into arteries and veins and finally AVR calculation. This chapter is mostly based in the paper "Assessment of Vascular Changes in Retinal Images" [12]. This chapter is organized as follows. In Section 5.1, the importance of the retinal vascular changes measurement, especially the AVR, is discussed. Section 5.2 presents the methods required for accomplishing the distinct phases of the AVR estimation process. The results of the tests on the images of INSPIRE-AVR database are presented in Section 5.3, where a comparison with the reference AVR values is also included. Finally, Section 5.4 summarizes the conclusions of this chapter. 85 86 Assessment of Retinal Vascular Changes 5.1 Background Retinal vessels are affected by several systemic diseases, namely diabetes, hypertension, and vascular disorders. In diabetic retinopathy, the blood vessels often show abnormalities at early stages [112], as well as vessel diameter alterations [6]. Changes in retinal blood vessels, such as significant dilatation and elongation of the main arteries, veins, and their branches [3,6], are also frequently associated with hypertension and other cardiovascular pathologies. Results of different studies indicate that retinal vasculature can be deeply affected by hypertension. There is a strong association between changes in retinal vascular and blood pressure in both adult and child populations [113,114]. Generalized retinal arteriolar narrowing is considered to be as an early characteristic sign of hypertension. The results of the Atherosclerosis Risk in Communities (ARIC) study in the United States showed that retinal arteriolar diameter is strongly and inversely related to higher blood pressure levels [115]. These data support the fact that there is a large connection between generalized arteriolar narrowing and hypertension. These studies have mainly focused on the use of smaller AVR as a way to measure generalized retinal arteriolar narrowing. Unfortunately, clinical examinations based on ophthalmoscopy are unlikely to be capable of detecting the subtle degree of arteriolar narrowing. Advances in medical image analysis, has made it possible to detect the smallest changes in arteriolar and venular diameters. Retinal image analysis can be used in order to study the cerebral microvasculature and related diseases due to the similarities of the retinal and cerebral vasculature in embryological origin, anatomical features, and physiological properties [116]. There is solid and consistent evidence that there is an association between retinal vascular changes and both clinical and subclinical stroke and a variety of cerebrovascular conditions independent of typical risk factors such as hypertension, diabetes, and smoking. The results of ARIC study showed that smaller retinal AVR values which is interpreted as generalized retinal arteriolar narrowing, can be an independent predictor of incident of stroke in middle-aged individuals. Researchers show that retinal microvascular flow is reduced in individuals with cerebral small vessel disease [116], and the retinal and cerebral arteriolar histopathology in these patients is similar to stroke cases. Discovery of the deficits related to the use of retinal AVR, subsequently encouraged newer studies to evaluate the association of retinal arteriolar and venular calibers separately with the risk of stroke [117,118]. By the use of this approach, two studies including the Rotterdam Eye Study [119] and the Cardiovascular Health Study [2] showed that, independent of 5.2 Arteriolar-to-Venular Ratio (AVR) Calculation Method 87 other stroke-related risk factors, wider retinal venular caliber predicted the future risk of clinical stroke rather than narrower retinal arteriolar caliber. As a result, the importance of arteriolar and venular caliber in prediction of stroke risk still remains to be determined. Using retinal vessel measurements, recent studies have tested the hypothesis that diabetes may have etiologic links with microvascular disease. The ARIC and Beaver Dam Eye studies [119,120], indicate that nondiabetic individuals with smaller AVR have higher risk of developing diabetes, independent of other diabetes-related risk factors. The results of Beaver Dam Study showed that this association is significantly stronger in individuals with hypertension at baseline [121]. The conductors of Rotterdam Eye Study suggested that this association may be a result of retinal venular dilatation rather than arteriolar narrowing. They verified an association of larger retinal venular caliber with pre-diabetes. This finding is consistent with earlier findings of prevalence data from the ARIC study [122]. While the responsible biological mechanisms for these observations remain to be clarified, experimental studies have shown that dilatation of retinal venules in normoglycemic patients can be achieved by administration of intravenous dextrose. Furthermore, reduced vascular reactivity associated with endothelial dysfunction and inflammatory processes may also play an essential role in the development of wider retinal venules and diabetes [123]. Several characteristic signs associated with vascular changes are measured aiming at assessing the stage and severity of some retinal conditions. Among the mentioned signs, generalized arteriolar narrowing, which is inversely related to higher blood pressure levels [117,118], is usually expressed by the Arteriolar-to-Venular diameter Ratio (AVR). The AVR value can also be an indicator of several diseases, like diabetes, hypertension and retinopathy of prematurity [124]. 5.2 Arteriolar-to-Venular Ratio (AVR) Calculation Method The estimation of AVR requires the detection of several retinal landmarks, namely the optic disc and the vessels, followed by accurate vessel caliber measurement and artery/vein classification [75,125]. Vessel segmentation is used for finding the vessels, and optic disc detection and segmentation is necessary to locate the ROI where vessel diameters are to be measured. An automatic AVR measurement system must also classify the retinal vessels into arteries and veins with high accuracy since small classification errors can have a significant influence on AVR values. Finally, caliber measurements are used for computing AVR according to the formula proposed 88 Assessment of Retinal Vascular Changes by Knudtson et al. [75]. Figure 5.1 shows the block diagram of the proposed method, which is detailed in the following subsections. Figure 5.1: Block diagram of the proposed method for AVR estimation. Vessel segmentation Vessel segmentation is used for finding the vessels, and as first stage for optic disc detection and A/V classification. Vessel segmentation is also necessary for the measurement of vessel diameters. For segmenting the vessels, the method previously proposed by Mendonça et al. [22] was chosen which was described in Section 3.1. The segmented vascular structure generated by this method for the high resolution image presented in Figure 5.2(a) is the binary image shown in Figure 5.2(b). Vessel caliber measurement Vessel calibers are the inside diameters of the blood vessels which are estimated using the distance transform of the binary vessel image. The result of this operation is a gray-level image in which the intensity value of each pixel, dpis the Euclidean distances from the considered pixel (p) and the 5.2 Arteriolar-to-Venular Ratio (AVR) Calculation Method 89 (a) (b) (c) (d) Figure 5.2: (a) Input image; (b) Binary vessel image result; (c) Distance transform result; (d) Centerline image. closest background pixel. Then the image is thinned in order to extract the centerline of vessels. For each vessel centerline pixel, the vessel caliber, vc(p), is simply estimated by doubling the distance values minus one (vc =2d−1). Afterwards, for the estimation of vessel caliber for each vessel segment, the average of all vessel calibers for its centerline pixels is calculated. The results of distance transform and centerline extraction are shown in Figure 5.2(c) and Figure 5.2(d), respectively. Region of interest determination The literature suggest the measurement of the AVR as the ratio between artery and vein widths measured in several circumferences centered at the optic disc. In most of the approaches, the region of interest (ROI) for the AVR calculation is considered as the ring area between the circumferences with diameter dOD and 1.5×dOD, where dOD is the optic disc diameter [66,75,125,126]. 90 Assessment of Retinal Vascular Changes Benavent et al. [127] considered the region between 0.5×dOD and 1.5×dOD and used three measurements in each vessel corresponding to the intersections between the vessels and the circles with diameters 0.5×dOD,dOD and 1.5×dOD. Pose et al. [128] took into account a region from dOD outwards and they considered several measurements for the same vessel over the region. In the manual method proposed by Hubbard et al. [125], the experts can move outside this region in order to measure the branches of an artery and they took into account only a measurement for each vessel found in the region between circles with diameters dOD and 1.5×dOD. Li et al. [126] proposed a tracking process to find all vessel segments in the ring area between to circles with dOD and 1.5×dOD diameters, while, Niemeijer et al. [66] calculated an AVR in six concentric regions around the OD and after that they computed the final AVR as the average of these values. In this work, the AVR is calculated from the calibers of vessels inside a ROI which is defined as the standard ring area within 0.5 to 1.5×dOD from the optic disc margin [75]. As a consequence, both the localization of the optic disc center (ODC) and its diameter are required for automating the AVR calculation. Two main options were considered for ROI definition: in the first approach, the ODC is estimated using an automatic methodology based on the entropy of vascular directions described in Mendonça et al. [103]. Then, the ODC is used as the center of the ROI, which is afterwards established considering a fixed disc diameter adapted to the size and field of view (FOV) of the image under analysis; the second approach delineates the OD border using SBF-based OD segmentation algorithm which was described in Chapter 3. Then for the ROI determination, the ODC and the disc diameter for each image are obtained by fitting a circle to the extracted OD boundary. The two approaches for ROI delimitation are illustrated in Figure 5.3. The initial estimate for the ODC location, obtained by the entropy-based method, is shown in Figure 5.3(a) and the ROI defined using a fixed OD radius is presented in Figure 5.3(c). Figure 5.3(d) shows the result of segmented OD boundary using SBF-based method. The approximation of the OD boundary by a circle is illustrated in Figure 5.3(e). Figure 5.3(f) refers to ROI definition after OD segmentation. Artery/vein classification In order to classify a vessel as artery or vein, two alternatives are used in this work: a supervised automatic graph-based A/V classification method and the unsupervised classification method described previously in Chapter 4. 5.2 Arteriolar-to-Venular Ratio (AVR) Calculation Method 91 (a) (b) (c) (d) (e) (f) Figure 5.3: (a) Retinal image with the ODC (black cross) detected using the method based on the entropy of vascular directions; (b) Circular OD boundary using a fixed OD radius of 180 pixels centered on the initial ODC in (a); (c) ROI for AVR calculation (delimited by the two green circles) with a fixed OD radius; (d) OD boundary using SBF-based method; (e) Approximation of the OD boundary by a circle (radius of 215 pixels); (f) ROI for AVR calculation (delimited by the two green circles) and the estimated optic disc margin (white circle). 92 Assessment of Retinal Vascular Changes The two approaches mainly differ in the methodology for assigning the final A/V class to each one of the labels resulting from graph analysis. In the supervised approach described in Section 4.1.4.1, the structural information provided by the graph is combined with the individual setting of A/V class for each vessel provided by a linear discriminant analysis (LDA) classifier and a set of intensity features extracted from the image. The result of A/V classification for the input image in Figure 5.2(a) is displayed in Figure 5.4(a), where the red color is used for representing arteries and veins are shown in blue. Figure 5.4(b) shows the arteries and veins found inside the ROI. The unsupervised approach described in Section 4.1.4.2, uses the red intensity of all vessel pixels and k-means clustering algorithm to assign the final A/V classes. The results of A/V classification for the whole image and for the main arteries and veins inside the ROI are shown in Figure 5.4(c) and Figure 5.4(d), respectively. AVR Computation Different approaches to estimate the AVR value have been proposed in the literature [125,129], where the AVR is calculated as the quotient between the estimated averages of artery and vein widths. Parr and Spears [129,130] and Hubbard et al. [125] obtained formulas which are widely used in medical studies [66,126,131]. In these methods, the AVR is computed as the quotient of the Central Retinal Artery Equivalent (CRAE) and the Central Retinal Venular Equivalent (CRVE) as follows AVR =CRAE CRVE (5.1) The CRAE and CRVE equivalents represent the relation among a vessel trunk and its two branches and they are computed iteratively using all vessels using Parr-Hubbard’s formulas [125, 129]. The Parr-Hubbard’s formulas for the calculation of CRAE and CRVE have been obtained theoretically and empirically and were afterwards revised by Knudtson et al. [75] . The revised formulas are independent of the number of selected vessels and only use the six main arteries and veins. In this work, the revised formulation is used for the calculation of CRAE and CRVE. These formulas correlate highly with those previously proposed by Parr-Hubbard formulas [125,129], with the advantages of being more robust against variability in the number of selected vessels and being independent of image scale. 5.2 Arteriolar-to-Venular Ratio (AVR) Calculation Method 93 (a) (b) (c) (d) Figure 5.4: (a) A/V classification result using the supervised graph-based method; (b) Main vessels inside the ROI (supervised AV classification and fixed OD radius); (c) A/V classification results using the unsupervised graph-based method; (d) Main vessels inside the ROI (unsupervised AV classification and OD segmentation). The CRAE and CRVE values based on Knudtson’s revised formulas are calculated as follows Arterioles :ˆ Wa=0.88∗wa12+wa221 2(5.2) Venules :ˆ Wv=0.95∗wv12+wv221 2(5.3) where wa1,wa2, and ˆ Waare, respectively, the width of the narrowest artery, the width of the widest artery, and the estimate of parent trunk for those arteries; wv1,wv2, and ˆ Wvhave similar meanings for veins. For computing the CRAE, the set with the six largest arteries inside the ROI is first selected. 100 Assessment of Retinal Vascular Changes (a) (b) (c) (d) (e) Figure 5.8: Scatter plots and regression lines between reference AVR values (a) between Observer 2 and reference; (b) between Niemeijer’s method and reference; (c) between Method 1 and reference; (d) between Method 2 and reference; (e) between between Method 3 and reference. 5.3 Results 101 (a) (b) Figure 5.9: (a) Number of subjects with matched classification between methods and reference; (b) Number of subjects with mismatched classification between methods and reference As it can be seen in Figure 5.9(a) and Figure 5.9(b), the number of matched and mismatched classifications for the proposed approaches are similar to the ones for the second observer in the case of considering different AVR values as a threshold. For example by consideration the AVR threshold value equal to 0.66, the number of matched classification between the reference and observer 2 is 26, while the number of matched classification between reference and method 1, method 2 and method 3 are 26, 26 and 24, respectively. 102 Assessment of Retinal Vascular Changes 5.4 Concluding Remarks In this chapter, three automatic approaches for calculating AVR in retinal images are described which are supported by a new approach for the optic disc segmentation and by new supervised and unsupervised techniques for the artery/vein classification. The proposed approaches were assessed in the 40 images of the INSPIRE-AVR dataset. For comparison purposes, first a fixed OD diameter is defined identical to the one established by other authors for this set of images. This solution has achieved a mean error of 0.05, identical to the one obtained by the second human observer. However, in order to increase the independence of AVR calculation from the particular characteristics of the images to be evaluated, two new solutions that are alternatives to the first one are described. One of the methods complements the OD detection algorithm with a segmentation approach that allows the estimation of the actual disc radius of the image under analysis, thus making the definition of the ROI for AVR calculation a fully automated procedure. The other method is an unsupervised classifier using intensity features whose results still needs to be combined with the structural information extracted from the graph representation for A/V classification purposes. This is also an important step towards automation because the classifier is naturally adapted to the intensity characteristics of each particular image. Although the unsupervised classification approach has produced a slight decrease in the global accuracy of vessel classification, this fact had no special impact on the result of AVR calculation because only the main vessel are normally involved in the calculation protocol, and no significant classification errors were observed for these vessels. The AVR values presented in Table 5.1 and Table 5.2 show that the proposed methods have a performance similar to those of human observers. The low errors and good correlation with reference AVR values are promising and demonstrate that described approaches have a high potential for clinical application. Chapter 6 Experimental Results This chapter presents RetinaCAD (Retinal Computer-Aided Diagnosis), an automatic system for analysing retinal images aiming at the estimation of features derived from retinal vasculature useful for the early detection and diagnosis of several systemic diseases. From the ophthalmologic point of view, this system is also relevant for assessing the implications of such diseases in the eye, as well as to evaluate the response of patients to specific therapeutic approaches. Several semi-automatic retinal image analysis systems have been proposed but the development of a fully automatic system for the assessment of vascular changes is still open. In this chapter, the image analysis tools available in RetinaCAD are presented, namely vessel segmentation, vessel width estimation, artery/vein (A/V) classification and optic disc segmentation. A pipeline of these methods, already described in the previous chapters, allows the computation of some relevant vessel related indexes, namely CRAE, CRAE and AVR, as well as various geometrical features associated with vessel bifurcations. The evaluation of the system on images of a new dataset from a local hospital shows a low failure rate and a significant correlation between the AVR calculated for distinct images of both eyes of the same patient. The obtained AVR values are also compared with the clinical information available for the subjects in this dataset. The images and clinical information in this study are anonymized and informed consent was obtained from all individuals. This chapter is organized as follows. Section 6.1 describes the RetinaCAD system and the tools available in the framework of the application. Experimental results using the images of a dataset from a local hospital are presented in Section 6.2. Finally, Section 6.3 summarizes the conclusions of this chapter. 103 104 Experimental Results 6.1 RetinaCAD System The RetinaCAD is a fully automatic system for the segmentation and classification of retinal structures and for the measurement of vascular features. This system can analyse optic disc centered retinal images with variable resolution and camera field of view. The implemented image analysis methods in this system are the ones that are introduced in the previous chapters. The developed system interacts with the user by means of a Graphical User Interface (GUI) which is shown in Figure 6.1. The implementation of the system was performed using MATLAB® R2013a version 8.1 (MathWorks Inc, Nattick, Massachusetts, USA). 6.1.1 Tools The main retinal image analysis tools provided by this application were detailed in the previous chapters. All tools are fully automatic and have shown a good performance when applied both individually and in combination. The list of available tools is the following: 1. Vessel segmentation 2. Vessel centerline extraction 3. Optic disc localization 4. Optic disc segmentation 5. Graph generation 6. Graph modification 7. A/V vessel classification 8. ROI determination The system has two options for using these tools: 1. Option to process each one of these tools as an individual task. 2. Single click option which uses these tools to compute some relevant vessel related indexes, namely CRAE, CRAE and AVR, as well as various geometrical features associated with vessel bifurcations. 6.1 RetinaCAD System 105 (a) (b) (c) Figure 6.1: RetinaCAD graphical user interface; (a) Main screen; (b) Second display; (c) Report. 106 Experimental Results 6.1.2 Measurements In addition to the measurement of vessel calibers, this application automatically obtains some retinal image characteristics and computes several vascular features, namely the CRAE, CRVE, AVR and bifurcation geometrical features. The list of measurements provided by RetinaCAD system is the following: 1. Diameter of the retinal image (in pixels) 2. Location of the optic disc center 3. Optic disc radius 4. Central Retinal Artery Equivalent (CRAE) 5. Central Retinal Venular Equivalent (CRVE) 6. Arteriolar-to-Venular Ratio (AVR) 7. Bifurcation geometrical features, namely bifurcation index, asymmetry ratio, diameter ratio, area ratio and Optimality. Computation of bifurcation geometrical features Alterations in the geometry of the retinal vascular network reflect a state of vascular dysfunction. Several features representative of the retinal bifurcation geometry have been evaluated in association with various systemic conditions [132,133]. The bifurcation geometrical features are obtained from the graph representation of the vascular network. In each bifurcation, we consider a trunk vessel and two branches, the largest and the smallest, with diameters d0,d1, and d2, respectively. The branching angles are the basic measurements related to a bifurcation. At each bifurcation, the branching angles and vessel segment calibers are used for deriving features which were defined before in [133], such as: • Bifurcation angle (θ): angle between the smallest and largest branches. • Bifurcation index (λ): quotient between the diameter of the smallest and largest branches. λ=d2/d1(6.1) • Asymmetry ratio (α): cross-sectional area of the smallest branch divided by that of the largest. α=d2 2/d2 1(6.2) 6.1 RetinaCAD System 107 • Diameter ratios (λ1,λ2): branches diameters divided by trunk diameter. λ1=d1/d0,λ2=d2/d0(6.3) • Area ratio (β): sum of the cross-sectional areas of two branches divided by that of the trunk. β=d2 1+d2 2/d2 0(6.4) • Optimality (ρ): deviation of the junction exponents from the optimal value of 3 based on branching physiological principles [134]. ρ=d3 0−d3 1+d3 21/3/d0(6.5) 6.1.3 Additional Features In order to increase interactivity and usability, additional features and tools are added to this application: 1. Single click option for the CRAE, CRVE and AVR estimation 2. Visualization of intermediate and final results 3. Ability to process an individual task 4. Display of vessel bifurcation/intersection (graph node) local and global features using a mouse click, such as: • Vessel segment caliber • Angle between vessel segments • Type of nodes (crossing, bifurcation, high curvature, meeting or terminal points) 5. Display of bifurcation geometrical features for each individual node 6. Availability of a follow mouse tool using a second display window for comparison purposes (Figure 6.1(b)) 7. Generation of reports in PDF or Excel file (Figure 6.1(c)) 8. Tools for manual modification of several calculated results, namely for changing • ODC location • OD radius value • ROI delimitation • Vessel segments A/V classification • Number of nodes in the graph 108 Experimental Results 6.2 Evaluation Results In this section, the experiments performed to test the RetinaCAD system for AVR computation are described and the main results are summarized. It is worth mentioning that the system requires no parameters setting for different datasets since all the parameters are adapted to the image resolution and camera FOV. The system is now evaluated on a new dataset from a local hospital which is described in the following subsection. 6.2.1 Material The system was evaluated using the images of a dataset from Centro Hospitalar São João (CHSJ dataset). This dataset contains 564 images from 141 subjects, where for each subject there are four images, two from the right eye and two from the left eye; the two images of each eye were acquired with two different FOV values (45° and 30°). All the images are optic disc-centered and the ones with 30° FOV have resolution of 1444×1444 pixels and the resolution of images with 45° FOV is 2196 ×1958 pixels. The average age of the subjects in this study is 60 in the range between 19 and 81 years. 53% of subjects are female while 47% are male. The clinical information is also available for this dataset, namely blood pressure, diabetic status, grades of diabetic, hypertensive and vascular retinopathy. The clinical information and images in this dataset are anonymized and there is informed consent from the individuals. The images were used for investigating the robustness of the system to the use of distinct images of the same subject, through the assessment of the correlation between measurements from images of the right and left eyes and from images of same eye with different FOV. Some of clinical information is also used for system validation. 6.2.2 Experimental Validation For the 564 images from CHSJ dataset that were analysed using RetinaCAD, the results for 27 images (less than 5% of all images) were not accepted as a consequence of errors in vessel segmentation, A/V classification and OD segmentation in images with no discernible OD or in images with severe pathological conditions. Figure 6.2 shows some images from this dataset where the RetinaCAD system fails to calculate the measurements correctly. Although some of these errors can be solved using the manual modification tool which is included in this application, all the images of the subjects where the automatic procedures have 6.2 Evaluation Results 109 (a) (b) (c) (d) (e) (f) Figure 6.2: Example images from CHSJ dataset where the RetinaCAD system fails. failed (11 subjects) were manually excluded. In the following, the results are obtained for 520 images from 130 subjects. Some results of A/V classification inside the ROI and obtained AVR values for images this dataset are shown in Figure 6.3. Although there are some differences between the results of A/V classification for the small vessels in the images of same subject, the obtained AVR values are similar. 116 Experimental Results 6.2.3 Clinical Validation The proposed AVR methodology included in RetinaCAD is compared with some clinical information. Table 6.5 shows the relation between the mean and standard deviation (SD) of AVR values for 130 subjects with available medical information. The AVR values for right and left eyes and also the average value of both eyes from the images with 30° FOV are compared with the medical status of subjects such as risk of diabetes and other pathological conditions. The subjects are categorized in different groups based on their pathological status. For each category, the mean and standard deviation of AVR values are obtained. The 95% confidence interval (CI) for the AVR values of each category are also included in table. The lower and upper limits of 95% CI of AVR values for each category are calculated using equations (6.6) and (6.7), respectively. Lower 95% limit =x−1.96×σ √n(6.6) Upper 95% limit =x+1.96×σ √n(6.7) where xis the mean AVR value, σrepresents the standard deviation and nis the number of subjects in each category. The AVR values for the categories of subjects with diabetes, diabetic retinopathy, hypertensive retinopathy and vascular retinopathy are compared in this table. In the last two rows of Table 6.5 the subjects are categorized in two groups based on all mentioned pathological conditions and their blood pressure (BP). The subjects with systolic BP higher than 140 mm Hg or diastolic BP higher than 90 mm Hg are considered as the ones with high blood pressure. The number of subjects with high BP is 41, while the rest of subjects have normal BP. Although there is a considerable overlap between the measurements, due to uncertainty of the measurement given by the SD, the average of AVR values for the categories of subjects with these diseases are smaller than the average of AVR values for the subject with non-pathological conditions. As can be observed in this table, the mean AVR value of 0.66 (95% CI: 0.64 - 0.68) in the nonpathological subjects with normal BP is higher than other group with pathological and/or high BP with mean AVR value of 0.62 (95% CI: 0.61 - 0.63), which shows evidence that the measurements 6.2 Evaluation Results 117 obtained for AVR are consistent with the expected association of the AVR value with the risk of systemic diseases and hypertension. Figure 6.8 shows the mean of AVR values and 95% confidence interval for each category of subjects. Red color represents the categories with pathological conditions and green color indicates the non-pathological categories. Figure 6.8(a) and Figure 6.8(b) represent the AVR values for right and left eyes, respectively, while the average of AVR values from both eyes for each category are illustrated in Figure 6.8(c). Table 6.5: Average of AVR values and 95% confidence intervals for the subjects with different pathological conditions. AVR (Mean ± SD) Number [95% CI] Pathological status of subjects Right eye Left eye Average of both eyes Diabetes 58 0.61 ± 0.08 0.62 ± 0.06 0.62 ± 0.06 [0.59 - 0.63] [0.60 - 0.63] [0.60 - 0.63] Diabetic retinopathy 21 0.63 ± 0.10 0.60 ± 0.06 0.62 ± 0.07 [059 - 0.67] [0.58 - 0.63] [0.59 - 0.65] Hypertensive retinopathy 63 0.61 ± 0.09 0.62 ± 0.07 0.61 ± 0.06 [0.59 - 0.63] [0.60 - 0.63] [0.60 - 0.63] Vascular retinopathy 69 0.61 ± 0.08 0.62 ± 0.07 0.61 ± 0.06 [0.59 - 0.63] [0.60 - 0.63] [0.60 - 0.63] Pathological or high BP 102 0.62 ± 0.08 0.62 ± 0.07 0.62 ± 0.06 [0.61 - 0.64] [0.61 - 0.64] [0.61 - 0.63] Non-pathological and normal BP 28 0.66 ± 0.06 0.65 ± 0.08 0.66 ± 0.06 [0.64 - 0.69] [0.62 - 0.68] [0.64 - 0.68] BP: Blood Pressure CI: Confidence Interval 118 Experimental Results (a) (b) (c) Figure 6.8: The average of AVR values with 95% confidence intervals for the subjects with different pathological status (a) Right eye; (b) Left eye; (c) Average of right and left eyes values. 6.3 Concluding Remarks 119 6.3 Concluding Remarks In this chapter, a user-friendly system, RetinaCAD is introduced which is capable of automatic detection, measurement and classification two main retinal landmarks, the optic disc and the vessels. RetinaCAD can measure several vascular features that are recognized as indicators for some prevalent systemic diseases. This application was assessed in the images of the CHSJ dataset, where it showed a relation between obtained AVR values and pathological status of 130 subjects. The lower AVR values in the people with diabetic, retinopathy conditions or high blood pressure in contrast to the nonpathological ones with normal blood pressure demonstrates the potential and usefulness of this system as a CAD tool for early detection and follow-up of diabetes, hypertension or cardiovascular pathologies. After comparing the measured AVR values for images of the same eye but with different FOV, a substantial correlation and a low mean error were achieved, thus allowing the conclusion that RetinaCAD is adequate both for research and for general clinical use. The evaluation of system on new datasets is still a challenge and there are still several issues and problems that need to be resolved. 120 Experimental Results Chapter 7 Conclusions and Future Work 7.1 Summary and Conclusions This thesis dealt with the development of automated retinal image analysis methods for the assessment of signs related with the changes in vessels calibers caused by several pathologies such as diabetes, hypertension, cerebovascular and cardiovascular diseases. Among several retinal vascular signs, Arteriolar-to-Venular Ratio (AVR) is a well known health biomarker and there is a strong need to develop an automated system for an accurate and reproducible estimation of AVR which requires different image analysis tools, namely vessel segmentation, vessel caliber estimation, optic disc (OD) segmentation and artery/vein (A/V) classification. This work was mostly focused on the development of methodologies for OD segmentation, A/V classification and automated calculation of the AVR value which are integrated in the RetinaCAD system with a potential for clinical applications. The starting point of this work was based on a vessel segmentation method previously developed by our research group. Some parts of this work are already published in peer-reviewed journals and presented in the different conferences. The list of related publications are included in Appendix B. AVR is estimated from the calibers of the vessels inside a specific region of interest (ROI), defined as the standard ring area around the OD. As a consequence, both the localization of the optic disc center and its diameter are required for automating the AVR calculation. For this reason, we proposed a new automatic method for OD segmentation using a multiresolution sliding band filter. The evaluation results on the images of three datasets demonstrate the better performance of the proposed method compared to recently published OD segmentation approaches and prove the independence of this method from changes in image characteristics such as size, quality and 121 122 Conclusions and Future Work camera field of view. Our approach is also able to produce useful results even in the presence of severe pathological conditions which shows the robustness of the SBF-based approach in the presence of diabetic retinopathy and risk of macular edema. We also presented a new automatic method for A/V classification based on the analysis of a graph extracted from the retinal vasculature. The proposed approach classifies the entire vascular tree through the combination of structural information with vessel intensity information, while most of the recent methods mainly use intensity features for discriminating between arteries and veins. Our method is evaluated using manual labeling for three public databases. The results demonstrate that this method outperforms recent approaches for A/V classification. On the other hand, the high accuracy achieved especially for the largest arteries and veins, confirm that this A/V classification methodology is reliable for the calculation of AVR and other indicators associated with vascular alterations. An automatic approach for the estimation of AVR in retinal images was introduced in this work. The proposed method includes new solutions for the OD segmentation and A/V classification. This method was evaluated using the images of the INSPIRE-AVR dataset. The low mean error of the measured AVR with respect to the reference ones were identical to the one achieved by a medical expert using a semi-automated system, thus demonstrating the reliability of the proposed solution for AVR estimation. Finally, the developed approaches are integrated in RetinaCAD, a system for the fast, reliable and automatic measurement of CRAE, CRVE, and AVR values, as well as several geometrical features of the retinal vasculature. RetinaCAD automatically identifies important landmarks in the retina, such as the blood vessels and optic disc, and performs artery/vein classification and vessel width measurement. The system is validated on a new dataset from local hospital and the results of AVR estimation shown a substantial correlation between images of same eye acquired with different camera fields of view. The clinical validation on the estimated AVR values showed a lower value in subjects with diabetes or pathological conditions in contrast with subjects with nonpathological conditions which gives some confidence on the potential of the system as a CAD tool for early detection and follow-up of diabetes, hypertension or other cardiovascular pathologies. However, we are aware that this needs extensive validation on much larger datasets in order to increase the degree of confidence and to have a good generalization of the results achieved so far. 7.2 Future Directions 123 7.2 Future Directions The overall investigation in this thesis has shown the usefulness of the RetinaCAD system and the proposed approaches for OD segmentation and A/V classification, as well as the methodology for the assessment of vascular changes particularly the computation of AVR values. However there is still great possibilities for the improvement of medical image analysis tools based on the needs and trends of different applications. Moreover, RetinaCAD still needs a more extensive validation on much larger datasets in order to increase the degree of confidence and to have a good generalization of the results achieved so far. As earlier described, a new graph-based method is proposed for the classification of vessels as arteries or veins. Since the extracted graph represents the vascular tree, it can be useful for identifying different vascular signs such as vascular bifurcation angles, branching patterns and fractal based features which can have significant impact on the early detection of several systemic diseases. On the other hand, with the development of image analysis methods for assessing retinal vascular changes, there is a large demand for fast and fully-automatic algorithms to register two or more retinal images which were taken at different times and from different views. Accurate and fast registration of retinal images is still a challenging problem since the low content contrast, large intensity variance as well as various pathologies caused deterioration in retinal images. For this reason, we expect that the graph representation of the vascular structure could be also a valuable contribution for the development of an efficient registration method. In the proposed unsupervised approach for the assignment of A/V classes to the labels of each subgraph, only one feature (red intensity) is used in the k-means clustering algorithm. It would be interesting to see the performance of this unsupervised technique using other features, such as green component, hue and intensity. From the clinical point of view, further studies are needed with more patients in a long duration follow-up diagnosis programme. In this work only the AVR values are used for the clinical validation and it would be interesting to investigate the relationship of obtained CRAE, CRVE and bifurcation geometrical features with clinical information. Furthermore, there are other retinal vascular signs which can provide clinically useful information to aid prevention and management of systemic diseases. Some of these retinal vascular signs are tortuosity, branching patterns, fractal geometrical features, focal arteriolar narrowing and arteriolar color changes. 124 Conclusions and Future Work We believe that RetinaCAD is a remarkable step toward developing an automated retinal computer-aided diagnosis system. One main goal as future work on the RetinaCAD system, can be the implementation of more retinal image analysis tools for the measurement of other retinal vascular signs which can provide several opportunities for further research. Appendix A Materials This Appendix summarizes the features of the retinal image datasets which are used for the evaluation of different retinal image analysis algorithms. These datasets were created by several research groups and are publicly available [105,110,111,137–141]. The specifications of the datasets that are described in this appendix are shown in Table A.1. Table A.1: Datasets specifications. Dataset name Number of images Image size FOV Available ground truth STARE 20 700 × 605 px 35° Vessel segmentation DRIVE 40 565 × 584 px 45° Vessel segmentation Artery/vein classification VICAVR 58 768 × 576 px - Artery/vein classification INSPIRE-AVR 40 2392 × 2048 px 30° AVR value Artery/vein classification Optic disc segmentation MESSIDOR 1200 1440 × 960 px 45° Risk of macular edema grade 2240 × 1488 px Diabetic retinopathy grade 2304 × 1536 px Optic disc segmentation ONHSD 99 760 × 570 px 45° Optic disc segmentation A.1 STARE Dataset The STARE dataset consists of 20 images for blood vessel segmentation [137]. These retinal images were captured using a TopCon TRV-50 fundus camera at 35° field of view, and afterwards digitized to 700×605 pixels, 8 bits per RGB channel. Two observers manually segmented all the images. 125