1 Introduction
Eye fundus retinal vessels are the only blood vessels of the body that can be visible with non-invasive imaging (Miri
et al.,
2017; Fraz
et al.,
2014). Retinal images are commonly used in manual and in automatic vessel structure analysis. This analysis leads to disease detection which is especially important in the early stages of the disease (Li and Wee,
2014). First clinical findings in early diabetic retinopathy are often made in abnormalities of retinal blood vessels (Miri
et al.,
2017). One of these abnormalities is the ratio between the artery and vein diameters (
AVR).
AVR is approximately equal to 0.7 for a healthy person (Sun
et al.,
2009). Ratio increases in the evaluation of sick persons, arteries become larger. Blood vessel width is measured for main vessels in an
AVR measurement region – the circular region from 1.5 to 3 optic nerve disc (
OD) radius from its centre (Knudtson
et al.,
2003). The ratio is calculated for vessels above and below the
OD centre separately. There are also some other methods for measuring this ratio, including more vessels for calculations, selecting other places for measurements, etc. (see e.g. Sun
et al.,
2009; Niemeijer
et al.,
2011). The timely evaluation of disorders can lead not only to the better evaluation of the disease but also to protection of disorder progression.
The measurement making and other eye fundus analysis steps are time consuming for an expert, and sometimes ophthalmologists cannot spend much time for analysis if they have many patients. Here, computer aided systems help to do most of the analysis and act as advisory or tracking systems (see, e.g. Li and Wee,
2014; Bankhhead
et al.,
2012). When making such systems, there are many specific situations to be considered. For example, when an expert does measurements, the measurement place is selected very differently for the specific image – the most important factor is the ability to make precise measurements. When the computer makes automatic measurements the place where they are done is always concrete, and in many cases, these places are not the best for measurements to make. The automatic system does not know the quality of the image and even whether the person is sick or not. In many cases (Muramatsu
et al.,
2011; Bankhhead
et al.,
2012) computer aided systems are adapted for specific technical equipment or even to the specific image database and when a new type of images is gathered, the system must update to meet the new requirements. Sometimes this leads even to an information loss because of reduced image resolution (Ravishankar
et al.,
2009).
Many eye fundus analysis algorithms are tested on high quality image databases only. However, in some specific circumstances, it is hard to get high quality images, or it is hard to obtain another image if the image is not of high enough quality. Image quality is affected by the illness also. These situations lead to the requirements for the automatic system to be able to deal with images of lower quality.
Classification is required for AVR calculations. In order to do this, main vessels must be detected and classified into arteries and veins. Then, with the view to evaluate the artery and vein ratio from eye fundus images at the specific place, the following steps can be applied. For example, OD should be detected in order to make measurements at a specific place. Then, places for measurements must be selected.
The
OD detection is a very important stage in eye fundus analysis (Buteikienė
et al.,
2012). As it was mentioned before,
OD detection is used for measurement place selection for
AVR calculations. It is also used in vessel tracing algorithms in disease evaluations based on
OD analysis like
OD excavation and in other cases. Main problems arising in precise
OD detection are uneven lightening, disease based changes in eye fundus and noise produced by vessels inside of
OD. In this paper, a two stage
OD detection algorithm is used. This
OD detection algorithm is based on the algorithm described in Stabingis
et al. (
2018) with several improvements described in the next section.
In the next stage, the vessel width measurement algorithm (Stabingis
et al.,
2018) is applied. Profile analysis is a technique used in vessel measurement (see e.g. Li
et al.,
2003; Bankhhead
et al.,
2012). The vessel width measurement algorithm, used in the proposed algorithm, is based on vessel profile analysis described in the paper (Stabingis
et al.,
2018). Profile information is often extracted according to the Bresenham algorithm which extracts information fast but makes distances between points distributed unevenly. In the proposed algorithm, profile information extraction uses spatial distance function which enables to use the same analysis algorithm with images of different resolution. Measurement algorithm is capable of detecting vessels without prior information about vessel existence at a specific point. The extracted vessel tree is used only for the detection of the initial measurement direction and only when this tree information is available.
AVR measurement zone vessels are measured, and then vessels classification is applied.
Classification problem is still very significant part in medical image analysis (see e.g. Miri
et al.,
2017; Renukalatha and Suresh,
2018; Morkunas
et al.,
2018). Many different classification algorithms are used for vessel classification in eye fundus images (see e.g. Miri
et al.,
2017; Dashtbozorg
et al.,
2014; Sun
et al.,
2009; Niemeijer
et al.,
2011; Treigys
et al.,
2008). If the classification is intended to be done automatically and with images of different type, unsupervised classification methods are more suitable. Supervised classification methods are more accurate than unsupervised methods, but they must be adapted for different type images. Different medical image analysis methods have specific image quality issues (see e.g. Prasath,
2017). Eye fundus vessel classification struggles with different or uneven lightening issues. The proposed algorithm uses
k-means clustering method with spatially recalculated features decreasing the influence of uneven lightening.
There are many different features which are commonly used in vessel classification. These features are often extracted from
RBG,
HSI, grey level image information using statistical functions, texture features and others (see e.g. Renukalatha and Suresh,
2018; Muramatsu
et al.,
2011; Stabingis
et al.,
2016). Such features are commonly used for many other image analysis problems. New type features, based on the profile information gathered from the vessel width measurement algorithm, are proposed in this paper. These features are extracted from the vessel measurement stage and use information from vessel inner part and average it along the vessel adaptively to the vessel width.
Classification accuracy is compared with other common classification algorithms of an eye fundus blood vessels according to publicly available eye fundus databases.
IVSIPIRE-AVR (Niemeijer
et al.,
2011) and
DRIVE (Staal
et al.,
2004) databases are used for comparison.
DRIVE database is used only for a demonstration that algorithm works with images of any resolution without the need for adaptation. The accuracy results from this database are not important, because such small resolution images are not used any more. The proposed algorithm is compared with algorithms proposed in papers (Dashtbozorg
et al.,
2014; Muramatsu
et al.,
2011; Mirsharif
et al.,
2013; Niemeijer
et al.,
2011). Algorithms used for vessel patch classification are selected only. Algorithms used for the whole vessel tree classification are not used for comparison, because our algorithm is built only for
AVR evaluation at this stage.
The method described in Muramatsu
et al. (
2011) is used especially for
AVR evaluation. This method consists of the following stages. First, vessel segmentation method is applied to the retinal image. The segmentation method involves top-hat transformation, double-ring filtering, thresholding techniques and it is performed on green colour channel. Vessel centreline pixels are extracted. According to vessel crossing and bifurcation points, centreline pixels are divided into separate vessel parts. Small vessels (with diameter smaller than 3 pixels) are removed from calculations. Thresholding and ellipse fitting is applied for
OD detection. The longer ellipse axis is used as
OD radius for
AVR measurement zone. Segmentation and
OD detection parameters were selected for specific image resolution. Centreline pixels from
AVR measurement zone were used in classification. For each centreline pixel, five features are extracted for classification. These features are
RGB values corresponding each pixel and contrast values gathered from red and green channels. The contrast features are calculated from red and green channels as the difference in average pixel values in a
$5\times 5$ pixel region around the pixel of interest inside the vessel, and in a
$10\times 10$ pixel region around the pixel of interest outside the vessel.
LDA classifier was used for each centreline pixel classification and each vessel segment is classified into an artery or a vein by selecting a majority class from classified centreline pixels. Then, major vessels are selected for
AVR evaluation. Vessel measurements are done according to segmented vessels, averaging along the vessel patches. This method was tested on
DRIVE image database and main results used for comparison are presented in Table
2.
The method used in Mirsharif
et al. (
2013) consists of the following steps: image enhancement, vessel segmentation, thinning, feature extraction, artery/vein classification and post processing. Histogram matching, several histogram equalization and multi scale retinex techniques were used for image enhancement. The method proposed in Soares
et al. (
2006) is used for vessel segmentation. This method uses Gabor wavelets. Vessels thinner than three pixels were removed from the vessel network. Vessels are divided into smaller segments and for each segment many features were extracted. Features included information from
RGB,
HSL and
LAB colour spaces and their statistics. From all these features, 8 best discriminating features were selected.
LDA classifier was used for this method. Then post processing, using structural knowledge, is applied for reclassification of the vessels unconnected to the tree. The method was used for classification of all vessels and for classification of vessels in
AVR measurement zone.
OD was marked manually here.
AVR measurement zone data from
DRIVE image database is used for comparison.
Method used in Niemeijer
et al. (
2011) consists of the following steps: pre-processing, vessel segmentation and centreline extraction,
OD detection, vessel width measurement, vessel classification,
AVR evaluation. Pre-processing involves
FOV extraction and Gaussian filtering. Vessel segmentation is performed with classification method. This stage was constructed for DRIVE database, and the use with larger images involves reducing the image size. Pixel classification produces likelihood map and threshold produces vessel structure. The method includes “tubogganing” technique from Fairfield (
1990). Then, the structure was thinned. A supervised position regression method was used to detect the centre of the
OD (Niemeijer
et al.,
2009). The
OD radius was considered as constant for all images. Vessel measurements were done perpendicular to detected local centreline angle. Artery/vein classification was done with
LDA method, using 27 different features. Features were extracted from different image channels using different statistics. Statistics were applied to vessel centreline pixels and to vessel patches. During this classification stage, vessel likelihood was calculated. Final class label for vessels was obtained only during the
AVR calculation stage. During this stage, vessels were classified into arteries and veins including prior artery/vein structural knowledge and classified vessel pairs were used for
AVR measurements. This method was tested with
INSPIRE-AVR database, which was first introduced in this paper (Niemeijer
et al.,
2011).
Finally, the method (Dashtbozorg
et al.,
2014), which includes graph analysis, is selected for the comparison. This method uses the characteristics of retinal vessel tree. The method described in Mendonca
et al. (
2013) is used for vessel segmentation. Vessel centreline structure is extracted and then used for graph construction. The obtained graph is then modified, to solve common error issues. Graph analysis leads to initial node classification. Then,
LDA classifier is applied with features extracted from different image channels and with their statistics. 19 features were used for this classification part. Whole vessel tree classification was performed there, also, calculation only for
AVR measurement region vessels are presented. This method was tested with
DRIVE and
INSPIRE-AVR databases.
The proposed algorithm is also tested on the high resolution image database OPTO-AVR described in Section 3. Here, images are acquired with Optomed OY digital mobile eye fundus camera Smartscope M5 PRO from soldiers before and after sports load. The database consists of 86 different images.
The structure of this article is as follows. Method implementation chapter describes the proposed method which consists of scale parameter evaluation, OD detection, vessel measuring, feature extraction, and classification. Method comparison with other similar methods is presented in the result section. Conclusions and further development opportunities are described in the conclusion section.
2 Method Implementation
Main adaptive features of the proposed method, compared with less adaptive methods are presented in Fig.
1. The proposed method is adaptive to different size images, no image downsampling is applied. Any parameters don’t need to be tuned manually for different image data sets, no different models, with different parameters are created. Image noise and uneven lightening of the image is considered while creating the algorithm. All main steps: OD detection, blood vessel measurements, artery/vein classification and AVR evaluation, are done automatically.
Fig. 1
Main adaptive features of the proposed method, compared with steps from less adaptive methods.
Main steps of an automatic algorithm for evaluation of artery and vein ratio are presented in Fig.
2. Green colour channel
${I_{G}}$ from eye fundus image is used for eye fundus mask
${I_{M}}$ extraction. For mask extraction, common image analysis techniques are applied. After thresholding, the higher intensity part from
${I_{G}}$, mathematical morphology closing and flood fill operations are used to form a circular mask of eye fundus part in the image.
${I_{M}}$ is then used to calculate the scale parameter
${p_{s}}$ according to Eq. (
1) combining logistic growth model with a linear model.
where
${\beta _{1}}$,
${\beta _{2}}$,
${\beta _{3}}$ are parameters of the logistic growth model part,
${\beta _{4}}$ is the parameter of the linear model part (
${\beta _{1}}=4$,
${\beta _{2}}=20$,
${\beta _{3}}=500$ and
${\beta _{4}}=3000$ in our case),
${W_{M}}$ is the width of fundus mask
${I_{M}}$. Parameter
${p_{s}}$, evaluated by Eq. (
1), is used in almost all steps of the algorithm. This parameter allows to use images of different sizes. Eq. (
1) is a modified parameter version from Stabingis
et al. (
2018) and includes adaptation for different
FOV values of images. Newer fundus cameras acquire images with larger resolution and also use optics with larger
FOV value. The obtained mask
${I_{M}}$ is further used in other steps of the algorithm like norming operation and others.
Fig. 2
Main steps of the proposed method for automatic AVR evaluation.
Image pre-processing allows to unify different images and to reveal some image features required for specific operations. Image pre-processing is treated as a vital stage for many eye fundus image analysis methods. The proposed method consists of several different large stages (see Fig.
2.), and for each stage, a different image pre-processing algorithm is applied. Pre-processing allows controlling different stages separately because some methods can enhance the extraction of the blood vessel tree but reduce
OD extraction reliability. Vessel tree
${I_{tr}}$ extraction stage is not very important in the proposed algorithm, because the
${I_{tr}}$ information is used in other stages only as the auxiliary. In some other eye fundus analysis algorithms,
${I_{tr}}$ extraction is very important, and the extracted tree is used as the reference for almost all other stages. Such dependence on extracted
${I_{tr}}$ can lead to some errors, especially for sick patients, when some artefacts in the image can be treated as blood vessels or when vessel and background contrast is very low, producing breakage in a vessel tree. There are several different technics for
${I_{tr}}$ extraction like thresholding, mathematical morphology, Gabor filtering, wavelets, etc. More precise methods are also developed. These methods produce higher extraction accuracy but require a larger amount of computations. Mathematical morphology based method is applied here and is described in more detail in Stabingis
et al. (
2018). All other stages work even without the extracted
${I_{tr}}$ and in places where
${I_{tr}}$ is incorrectly extracted.
2.1 Optic Nerve Disc Extraction
OD extraction algorithm consists of two stages. Initial
${\mathit{OD}_{\mathit{init}}}$ is located at the first stage, and then real
${\mathit{OD}_{r}}$ circle is detected and evaluated. During the second stage, the circle is fitted to detect a circle like an object in the image. Full
OD detection scheme is presented in Fig.
3. Main
OD detection steps are described in more detail in Stabingis
et al. (
2018).
Fig. 3
Two stage OD detection scheme. Green dashed lines show the influence of scale parameter, used for image size adaptation.
First
${\mathit{OD}_{\mathit{init}}}$ detection stage detects the preliminary location of
OD. During this stage the probability map
${I_{\mathit{pr}}}$ is created. This probability map combines the influence of five different probabilities: image information part
${I_{\mathit{pr}}}$, vessel tree line intersection map (Stabingis
et al.,
2018) multiplied by 0.2, vertical centre probability map multiplied by 0.2, vertical gradient map multiplied by −0.4, and the horizontal gradient map multiplied by −0.3.
${I_{\mathit{pr}}}$ is used for finding the most intense circle
${\mathit{OD}_{\mathit{init}}}$ with a radius from
$80\cdot {p_{s}}$ to
$200\cdot {p_{s}}$ intersecting with
${I_{tr}}$.
Normally eye fundus image is photographed in such a way that
OD is located near the vertical centre of the image. In order to add this feature, vertical centre probability map is calculated according to Eq. (
2).
where
h is
${I_{I}}$ image height,
${\sigma ^{2}}={(\frac{h/2}{3})^{2}}$,
${\beta _{v}}$ is intensity control parameter for vertical probability map (
${\beta _{v}}=100$ in our case),
i is image line number and
j is image column number. To eliminate uneven lightening vertical and horizontal gradient maps are calculated. A vertical gradient map is calculated by averaging every image row, and a horizontal gradient map is calculated by averaging every image column.
The region of interest is selected according to
${\mathit{OD}_{\mathit{init}}}$ circle multiplied by 2.5, and it is used for real
${\mathit{OD}_{r}}$ detection (Stabingis
et al.,
2018). After several image analysis and processing operations (see Fig.
3) Hough circle detection is applied, and ten best circles are selected. A circle surrounded by most points from
${I_{b}}$ is selected as
${\mathit{OD}_{r}}$ and is used in further analysis.
2.2 Blood Vessel Measurements
After
${\mathit{OD}_{r}}$ is detected, the information is used for vessel measurements. Measurements are done along the circle with radius
$2\cdot {r_{\mathit{OD}}}$ from
${\mathit{OD}_{r}}$ centre (
${x_{\mathit{OD}}}$,
${y_{\mathit{OD}}}$), where
${r_{\mathit{OD}}}$ is
${\mathit{OD}_{r}}$ radius. Measurements are done at five pixel length intervals calculated according to Eq. (
3).
where
${\alpha _{i}}={\alpha _{i-1}}+\frac{5}{2\cdot {r_{\mathit{OD}}}}$,
$alph{a_{i}}=0$ and
${\alpha _{1}}<2\pi $. Algorithm proposed in Stabingis
et al. (
2018) is used for vessel width measurements. The algorithm automatically detects the place where the vessel can be. After the vessel is detected and measured the algorithm goes along the detected vessel and measures other profile widths. Vessel data is collected for every vessel in the range of radiuses from
$1.5\cdot {r_{\mathit{OD}}}$ to
$3\cdot {r_{\mathit{OD}}}$ (see Fig.
5). Average vessel width value is used for further calculations. Measurements are omitted at places where the vessel is already measured. Vessel measurements are done by analysing vessel profile data. Profile information is gathered from the image using spatial distance based function Eq. (
4). For every profile, 100 intensity values are calculated from
J nearest pixels.
where
${v_{j}}$ is the intensity value for step
j at location
${s_{j}}$,
${h_{ij}}=1-\frac{d({s_{i}},{s_{j}})}{3}$ is the distance function, for
$d({s_{i}},{s_{j}})<3$,
$d({s_{i}},{s_{j}})$ is the Euclidean distance between
${s_{i}}$ and
${s_{j}}$ points. Profile information is smoothed with Gaussian filter and analysed for width evaluation. Profile analysis points are shown in Fig.
4. Profile minimum point is found near the profile centre. Then, another close minimum point is founded in order to eliminate possible central reflex information. Two points
${a_{l}}$ and
${a_{r}}$ are used as left and right profile parts. Local maximum points
${b_{l}}$ and
${b_{r}}$ are found. Then fastest decaying points
${c_{l}}$ and
${c_{r}}$ are detected. Similarly to
${c_{l}}$ and
${c_{r}}$, decaying marginal points
d are selected for left and right sides. Middle points
${e_{l}}$ (between
${d_{lb}}$ and
${d_{lt}}$) and
${e_{r}}$ (between
${d_{rb}}$ and
${d_{rt}}$) are selected as vessel width measurement points.
Fig. 4
Blood vessel profile analysis.
Such profile analysis is performed for profiles extracted for point (
${x_{i}}$,
${y_{i}}$) at different angles and the angle with the smallest profile width is selected as profile perpendicular to vessel. For the first measurement of the vessel, angles from
${0^{\circ }}$ to
${180^{\circ }}$ are analysed, and the detected angle is used for further vessel measurements decreasing the amount of different angles to be analysed. Vessel measurement algorithm measures vessels similar to the expert measurements and is presented in more detail in Stabingis
et al. (
2018).
2.3 Feature Extraction and Classification
Eye fundus vessel classification is a complicated task. For the proposed algorithm, novel features are extracted for classification. After analysing many different vessels from different type images it was considered that the best part for discriminating veins from arteries is the vessels inner part. Usually veins are darker than the arteries and arteries often have a brighter centre reflex line. The outer part of different type vessels can be very similar. Three different features are extracted from vessels for classification. These features are an average value between
${a_{l}}$,
${a_{r}}$ points, an average value between
${d_{lb}}$,
${d_{rb}}$ points and an average value between
${d_{lb}}$,
${a_{l}}$, and
${a_{r}}$,
${d_{rb}}$ points (see Fig.
4). Features are calculated from green image channel
${I_{G}}$ and from red channel
${I_{R}}$ separately. These two channels are used for their best discrimination information. Average feature value is calculated along the vessel part. So, for every detected vessel, six different averaged features are calculated.
As it was mentioned in the introduction chapter, uneven image lightening can lead to misclassification of vessels. In order to eliminate this influence, spatial distance based normalization function Eq. (
5) is applied to extracted features.
where
${\hat{\mu }_{d}}$ is the local mean and
${\hat{\sigma }_{d}^{2}}$ is the local variance of corresponding vessel features,
ω is the set of 10 closest vessels,
${d_{i}}$ is Euclidean distance between vessel patch centres.
Fig. 5
Result image from
INSPIRE-AVR database (Niemeijer
et al.,
2011) showing main analysis elements. Blue colour – classified veins, green colour – classified arteries. Main vessels used for
AVR evaluation are marked.
After distance based normalization, features are used for vessel classification. The
k-means clustering method is used. Top vessels and bottom vessels according to
${\mathit{OD}_{r}}$ centre point are classified separately. Vessel parts used for features are brighter for arteries than for veins, so a class with a larger feature values is treated as an artery and a class with smaller feature values as vein classes. Two largest vessels from different classes are selected for the top part and for the bottom part. According to the selected vessel widths, top and bottom
AVR values are evaluated. One eye fundus image with classification results and main analysis elements is presented in Fig.
5. Classification results and comparison are presented in the results section.