Int J CARS DOI 10.1007/s11548-013-0965-9
ORIGINAL ARTICLE
Delineation of blood vessels in pediatric retinal images using decision trees-based ensemble classification Muhammad Moazam Fraz · Alicja R. Rudnicka · Christopher G. Owen · Sarah A. Barman
Received: 15 July 2013 / Accepted: 18 November 2013 © CARS 2013
Abstract Purpose Automatic segmentation of the retinal vasculature is a first step in computer-assisted diagnosis and treatment planning. The extraction of retinal vessels in pediatric retinal images is challenging because of comparatively wide arterioles with a light streak running longitudinally along the vessel’s center, the central vessel reflex. A new method for automatic segmentation was developed and tested. Method A supervised method for retinal vessel segmentation in the images of multi-ethnic school children was developed based on ensemble classifier of bootstrapped decision trees. A collection of dual Gaussian, second derivative of Gaussian and Gabor filters, along with the generalized multiscale line strength measure and morphological transformation is used to generate the feature vector. The feature vector encodes information to handle the normal vessels as well as the vessels with the central reflex. The methodology is evaluated on CHASE_DB1, a relatively new public retinal image database of multi-ethnic school children, which is a subset of retinal images from the Child Heart and Health Study in England (CHASE) dataset. Results The segmented retinal images from the CHASE_ DB1 database produced best case accuracy, sensitivity and specificity of 0.96, 0.74 and 0.98, respectively, and worst case measures of 0.94, 0.67 and 0.98, respectively. Conclusion A new retinal blood vessel segmentation algorithm was developed and tested with a shared database. The M. M. Fraz (B) · S. A. Barman School of Computing and Information Systems, Faculty of Science Engineering and Computing, Kingston University London, London, UK e-mail: [email protected] A. R. Rudnicka · C. G. Owen Division of Population Health Sciences and Education, St. George’s, University of London, London, UK
observed accuracy, speed, robustness and simplicity suggest that the algorithm may be a suitable tool for automated retinal image analysis in large population-based studies.
Keywords Medical imaging · Retinal image analysis · Ensemble classification · Image segmentation · Boot strapped decision tree
Introduction The human retina has the potential to reveal important information about retinal, ophthalmic, and even systemic diseases such as diabetes, hypertension, and arteriosclerosis [1]. Retinal imaging is increasingly used in large scale, populationbased studies for detection of glaucoma, diabetic retinopathy, age-related macular degeneration and cardiovascular disease [2]. A number of retinal blood vessel features such as arteriolar nicking and narrowing have been linked to systemic disease, and the morphological characteristics of retinal blood vessels themselves have been discussed in association with cardiovascular and coronary disease in adult life [3] and with retinopathy of prematurity in infancy [4]. The morphology of retinal vessels (particularly arterioles) has also been examined in relation to cardiovascular risk factors both in early and adult life [4]. Retinal vessels often have strong light reflexes along their centerline, which is more apparent on arteriolars than venules, and in younger compared with older participants, especially those with hypertension. The morphological characteristics of retinal images of premature infants and children are very different than that of the adult retina. Choroidal vessels are more visible alongside the retinal vessels in retinal images taken from premature infants [5]. A bright central reflex on the vessels and illumination artifacts
123
Int J CARS
Fig. 1 Vessel profiles: a–c Vessel profile of image from DRIVE dataset d–f Vessel profile of image from CHASE dataset: First column is the retinal image, second column is the magnified view of image part
marked in first column, and third column is the mesh plot of intensity profiles of second column
contribute to challenges in image processing when retinal images from school children are considered [4]. A variety of blood vessel segmentation algorithms is present in the literature [6]. However, most of them are evaluated on a limited range of retinal images of adult people, which includes twenty images each from the DRIVE [7] and the STARE [8] databases, which do not cater for the image-related characteristics such as, inter-image and intra-image variability in luminance, drift in contrast and uneven background gray level values, which are often found in the images acquired in larger population-based studies. The vessel profiles of a normal healthy adult retinal image and the image of a child retina from the CHASE [9] dataset are illustrated in Fig. 1, where the low contrast of blood vessels and the uneven background with augmented noise is clearly visible. Detecting blood vessels in retinal images of school children with the presence of uneven illumination, non-homogenized background, bright and dark lesions and central vessel reflexes is a challenging problem. This paper presents a supervised method for segmentation of blood vessels in retinal images of 9- and 10-yearold children of different ethnic origin using an ensemble classifier of bagged decision trees. The images named as CHASE_DB1 [9] are characterized with stark differences in background levels of retinal pigmentation (being more pigmented in South Asians and African Caribbean compared
with white Europeans). The classifier utilizes a feature vector based on the vessel-ness measure obtained by utilizing a group of the dual Gaussian matched filter, the Gabor filter and second-order derivative of Gaussian matched filter along with the multiscale line strength measure of blood vessels. The feature vector encodes information to successfully handle the blood vessels containing the central vessel reflex and without it. The classifier based on the boot strapped decision trees is a classic ensemble classifier which has been broadly applied in many application areas of image analysis, but has not been extensively utilized for retinal image analysis aimed at retinal vessel segmentation. The estimation of classification accuracy and prediction of feature importance during the training phase without using test data are the significant features of the boot strapped ensemble classifier. Moreover, the algorithm is computationally fast in training in addition to classification because it needs fewer samples for training as compared to other supervised algorithms available in the literature. The organization of the paper is as follows. In Sect. 2, the retinal vessel segmentation methodologies available in the literature are reviewed. The methodology, the feature vector and the classifier are explained in Sect. 3. In Sect. 4, the performance metrics, the accuracy and the robustness of the algorithm is presented and discussed. Finally, the paper is concluded in Sect. 5.
123
Int J CARS
Related work There is a substantial amount of work reported in the literature on automated segmentation of blood vessels in retinal images. A recent detailed review of these methods can be found in [6]. The vessel segmentation algorithms can be divided into two broad categories: rule-based methods and supervised methods. The rule-based methods can further be classified in to the techniques based on matched filtering, morphological processing, vessel tracking, multiscale analysis and model-based algorithms. The matched filtering methodology employs the piecewise linear approximation, the decrease in vessel diameter along vascular length and the Gaussian like intensity profile of retinal blood vessels and uses a kernel based on a Gaussian or its derivatives to enhance the vessel features in the retinal image. The matched filter was first proposed by Chaudhuri [10] and later adapted and extended by Hoover [11], Jiang [12] and Al-Rawi [13]. Mathematical morphology in combination with curvature evaluation [14] and matched filtering for centerline detection [15–17] is also exploited for retinal vessel segmentation. The algorithms based on vessel tracking [18,19] segment a vessel between two points using local information and work at the level of a single vessel rather than the entire vasculature. The multiscale approaches for vessel segmentation are based on scale-space analysis. The width of a vessel decreases gradually across its length as it travels radially outward from the optic disc. Therefore, the idea behind scale-space representation for vascular extraction is to separate out information related to the blood vessel having varying width. The multiscale second-order local structure of an image (Hessian) is examined and a vessel-ness measure is obtained on the basis of eigenvalue analysis of the Hessian [20]. A likelihood ratio test for vessel centerlines extraction that combines matched filter responses, confidence measures and vessel boundary measures was also presented [21]. The vessel cross-sectional intensity profiles can be approximated by a Gaussian or its higher-order derivatives, a mixture of Gaussians or hermite polynomial profile in cases where the central vessel reflex is present. Vermeer [22] modeled the vessel profile as a Laplacian kernel for segmentation. Mahadevan [23] presented a robust and modular framework for vessel detection in noisy images which is based on the estimation of the log likelihood of vessel parameters in a noisy environment. The framework is adaptable to incorporate a variety of vessel profile models including Gaussian, derivatives of Gaussian and dual Gaussian and various noise models like Gaussian noise and Poisson noise. A multiresolution hermite model was proposed for vascular segmentation by Li [24], which employed a two-dimensional hermite function intensity model in a quad-tree structure over a range of spatial resolutions. Lam et al. [25] proposed a vessel segmentation algorithm for pathological retinal images based on
the divergence of vector fields. The algorithm presented by Lam [26] was based on regularization- based multiconcavity modeling and was able to handle both normal and pathological retinas with bright and dark lesions simultaneously. Paulus et al. [27] proposed a methodology for quality assessment of retinal photographs. The cross-sectional intensity profile of retinal vessels was estimated by the dual Gaussian model initially used by Gao [28,29]. A number of authors have investigated the use of active contour models, or snakes [30] and geometric models based on level sets [31]. The supervised methods utilize ground truth data for the classification of vessels based on given features. These methods include the use of neural networks [32,33], principal component analysis [34], k-nearest neighbor classifiers [35], ridge-based primitives classification [36], Gaussian Mixture Model [37], line operators and support vector machine (SVM) classification [38], feature-based Ada-Boost classifier [39], SVM-based semi-supervised classification [40]. A decision tree-based ensemble classifier is presented in [41] for vessel segmentation which works well on normal and pathological retinal images. This work utilizes the ensemble classification approach aimed at segmenting the retinal vasculature from retinal images acquired from children attributed with normal vessels as well as the vessels with a central reflex.
The methodology The feature vector is designed to compute the quantifiable measure of the vessel-ness for each pixel in such a way that the ensemble classifier of bagged decision trees can successfully differentiate between the vessel (with and without the central vessel reflex) and the non-vessel region, an extension of work presented in [42]. A decision tree-based ensemble classifier is proposed in [41], where the main goal is to segment the blood vessels from adult retinal images containing pathologies. The color analysis of anatomical structures in retinal images reveals that the retinal pathologies appear either as dark lesions (e.g., hemorrhages) or bright lesions (soft and hard exudates). The utilized feature vector encodes information to successfully handle both normal and pathological retinas with bright and dark lesions simultaneously. In the current paper, we aim to solve the problem of retinal vessel segmentation from pediatric retinal images which are characterized with a bright central reflex in arterioles and wider venules; low contrast and uneven illumination. We have employed the feature vector which successfully differentiates the blood vessels containing the central vessel reflex and without it from the retinal background exhibiting poor contrast and uneven illumination. For this purpose the following changes are introduced in the feature vector.
123
Int J CARS
• A generalized multiscale line detector is used to compute the measure of the vessel-ness for each pixel in the retinal image. Besides highlighting the vessel structure, the generalized multiscale line detector reduces the false positive responses given to closely located vessels and at vessel crossovers, it deals with the central vessel reflex in an effective way. • The response of the Gaussian matched filter to the blood vessels with a central vessel reflex appears incorrectly as two blood vessels running parallel to each other. A twinGaussian model is used to give a correct response and imitate the vessels with a central light reflex which is apparent in the retinal images of children in the CHASE_DB1 database. • The retinal images in CHASE_DB1 contain a variation in background intensity due to non-uniform illumination. Therefore, we aim to obtain a shade corrected image by estimating and removing the background lighting variations using a simple linear transformation. The pixel intensity of the shade corrected image is included as one feature in the feature vector.
Moreover, we have also included a linear combination of line operator responses obtained at multiple scales, the morphological top-hat operation and gradient orientation analysis of the gradient vector field. The feature space is normalized to zero mean and unit standard deviation by applying a normal transformation, i.e., the filter response at each scale is subtracted by its mean and divided by its standard deviation. The inherent strength of ensemble classification for vessel segmentation is one in which multiple classifiers utilizing multiple cues/features for classification are strategically generated and combined to improve the classification or prediction performance of a model, or reduce the likelihood of a poor or unfortunate selection. This approach is intuitively used in our daily lives where we seek the guidance of multiple experts, weigh and combine their views in order to make a more informed and optimized decision. The boot strapped ensemble can estimate the classification accuracy and the feature importance index during the training phase without using the test data. The set of filters contain multiscale filters; therefore, the importance of the filter response at each scale is computed during classifier training, and the most significant features are selected to be included in the feature vector, thus reducing the dimensionality which in turn renders better computational performance. Moreover, the optimal number of decision trees used in ensemble creation can also be estimated during the training phase. The filters set and the line operator are applied to the inverted green plane of the colored RGB retinal image without pre-processing because the green channel of the RGB color fundus images gives the highest contrast between vessel and background. The green plane is
123
normalized and shade corrected in order to remove the variation in background illumination and shade correction before applying the morphological transformation and computing the Gradient Orientation Analysis (GOA) map. The feature vector Dual Gaussian and second derivative of Gaussian matched filter The small curvatures of blood vessels can be approximated by piecewise linear segments. The vessels appear in contrast relative to the background, have lower reflectance compared with the non-vessel retinal surface and almost never have ideal step edges. The intensity profile of the vessels can be approximated by a Gaussian curve and varies slightly from vessel to vessel. More specifically, the first-order derivative of Gaussian matched filter acts as the edge filter of the retinal vessel wall [6] and the second-order derivative of the Gaussian matched filter is employed for enhancement of piecewise linear structure of retinal vessel [43]. The Gaussian shaped matched filter and its extension have been used for vessel detection in retinal images [6]. The profile of the filter which is designed to match that of a blood vessel is defined as 1 x 2 exp − 2 ; for |x| f (x, y) = h 1 √ 2σ1 2π σ1 ≤ (t.σ1 ), |y| ≤ L/2
(1)
where, σ represents the scale of the filter; L” is the length of the neighborhood along the y-axis; h 1 is constant representing the height of the Gaussian curve, x = x cos θ + y sin θ, θ is the rotation angle of the filter kernel and t is a constant and is usually set to 3 because more than 99 % of the area under the Gaussian curve lies within the range [−3σ , 3σ ]. The response of the Gaussian matched filter to the blood vessels with a central vessel reflex looks like two blood vessels running parallel to each other. A twin-Gaussian model is used to imitate the vessels with a central light reflex which is apparent in the retinal images of children in the CHASE_DB1 database. For this purpose, a second small Gaussian curve is subtracted from the main Gaussian curve. The second Gaussian is oriented at the same angle and controlled by same parameters as the first one. The parameters of the second Gaussian control the height and width of central light reflex and the matched filter kernel is expressed as, 1 x 2 exp − 2 f (x, y) = h 1 √ 2σ1 2π σ1 1 x 2 exp − 2 − h2 √ (2) 2σ2 2π σ2
Int J CARS
Gaussian Matched Filter Kernel
Dual Gaussian Matched Filter Kernel
h1
h1
-X
0
+X
(a)
(b) Fig. 2 Dual Gaussian matched filter; a, b cross section and 2D profile of Gaussian matched filter, c, d cross section and 2D profile of dual Gaussian matched filter. (For the purpose of demonstration, the filter
-X
0
+X
(c)
(d) are constructed with following parameters; h 1 = 1, h 2 = 0.4∗ h 1 , σ1 = 2.5, σ2 = 0.8∗ σ1 , L = 27, t = 3 and = 0 degrees)
Fig. 3 Filter set responses; a retinal sub-image, b Dual Gaussian filter response highlighting the central reflex pixels, c second-order derivative of Gaussian matched filter response
In Fig. 2, the first column in the figure shows the cross section and two-dimensional plots of the Gaussian Matched filter kernel whereas the second column illustrates the cross section and two-dimensional profile of the dual Gaussian matched filter kernel. For the purpose of demonstration, the filter kernels shown in Fig. 2 are constructed with following parameters; h 1 = 1, h 2 = 0.4∗ h 1 , σ1 = 2.5, σ2 = 0.8∗ σ1 , L = 27, t = 3 and = 0 degrees. The second-order derivative of the Gaussian matched filter has also been exploited for vessel enhancement [43]. We have constructed a set of twenty four kernels in total,
consisting of two sets of twelve 2D dual Gaussian kernels and twelve second-order derivative of Gaussian kernels. The kernels in each set are oriented at twelve equiangular rotations θ1 · · · θ12 . The orientation angles of the kernels in each set of the filter group can be defined as {θ |0 ≤ θ ≤ π & θ mod(π/8) = 0}. The filter set is applied to the retinal image and the response is taken as the value for the highest scoring kernel at each pixel, which is taken as the pixel feature vector. Figure 3 illustrates a section of retinal image, the response of the dual Gaussian filter kernel and the secondorder Gaussian kernel. It is observed that the central vessel
123
Int J CARS
The retinal vasculature is composed of arteries and veins appearing as piecewise linear features, with variation in width and their tributaries visible within the retinal image. The concept of employing line operators for detection of linear structures in medical images is introduced in [44] which is modified and extended in [38,41,45] to incorporate the morphological attributes of retinal blood vessels. The average pixel intensity is evaluated along lines of a particular length passing through the pixel under consideration at 12 different orientations spaced by 15 degrees each. The line with the highest average pixel intensity is selected. The line strength of a pixel is calculated by computing the difference between the average gray values of a square sub-window centered at the target pixel with the average intensity of the selected line. A square sub-window of size w is centered at each pixel w . and average of pixel intensities is calculated, termed as Savg The average pixel intensity is measured along lines of a pixel length w passing through the center at 12 different orientations from 0 to 180 degrees at angular resolution of 15 degrees. The line with the highest average pixel intensity is w . The measure of the candidate line and is termed as Dmax vessel-ness of a pixel Vw is calculated by computing the difference between the average pixel intensity of a square subwindow centered at the target pixel with the average intensity of the candidate line as,
detector, the window size w is chosen in such a way that if the square sub-window is placed on the image and the center of this window corresponds to the center of vessel, then this window contains approximately the same number of vessel and background pixels. Therefore, it is selected to be twice the size of a typical vessel width in the image database. It has been experimentally found that W = 25 is a good choice for CHASE_DB1. The intensity of pixels in the middle of a vessel with a central reflex is usually lower than their surrounding pixels, instead of attaining the maximum value as in case of normal vessel. The line detector with longer lines can successfully recognize them as a vessel pixel because the candidate line includes only a small number of central reflex pixels; hence, the average intensity of pixels in the candidate line is not substantially affected and the central reflex pixels have a high vessel-ness measure. The vessel-ness image is shown in Fig. 4. It can be observed in Fig. 4c that the line detector with longer lines, i.e., a larger window size is effective in dealing with vessels with a central reflex but at the same time it produces false vessel responses at background pixels near the strongly contrasted vessels, i.e., those vessels which appear brighter than background. Also, it tends to merge the vessels in close proximity to each other and produces false responses in the area between vessels that are close to each other. These drawbacks can be avoided by using a generalized multi-scale line detector, which uses variable length of aligned lines in a fixed square sub-window. The vessel-ness measure using the generalized multi-scale line detector is can be expressed as;
w w − Savg Vw = Dmax
L w − Savg VwL = Dmax
reflex pixels are more highlighted with the dual Gaussian filter. Line strength features
(3)
The principal is that if the target pixel belongs to the vessel, then the vessel-ness measure will be large due to alignment of the candidate line along the direction of vessel. Contrary to that, the vessel-ness measure is low for a background pixel because the difference between the pixel intensity of the candidate line and square sub-window is small. In this basic line
(4)
where, 1 ≤ L ≤ w. The multi-scale line detector can be obtained by using different values of L. Besides highlighting the vessel structure, the shorter length line detector will not give false responses to closely located vessels and at vessel crossovers. Whereas, as described above, the longer length line detectors can deal with the central vessel reflex in an
Fig. 4 Line strength feature; a segment of retinal image; b, c Vessel-ness image with W = 3 and W = 21; d linear combination of multiscale line operator responses
123
Int J CARS
Fig. 5 Feature vector; a retinal image, b line strength feature, c Gabor filter response, d second derivative filter response, e gradient orientation analysis map, f morphological transformation
effective way. With the purpose of preserving the strength and exterminating the shortcoming of each individual line detector, the final vessel-ness measure is obtained by a linear combination of line responses computed by using line detectors at different scales. The combined vessel-ness measure can be defined as; 1 L Vw + I (5) Vcombined = nL + 1 L
where, n L is the number of scales used, VwL is the vessel-ness measure obtained using a line detector with line length L and square sub-window W , and I is the inverted green channel image. The combined vessel-ness image is shown in Figs. 4d and 5b, where it can be observed that a linear combination of responses from multi-scale line detectors works effectively in removal of background noise, dealing with the central vessel reflex while avoiding merging of close vessels and false responses near vessel crossovers. Multi-scale Gabor filter A Gabor filter is a linear filter and has been broadly used for multi-scale and multi-directional edge detection. The Gabor
filter can be fine-tuned to particular frequencies, scales and directions and therefore acts as low level feature extractor and background noise suppresser. The impulse response of a Gabor filter kernel is defined by the product of a Gaussian kernel and a complex sinusoid. It can be expressed as
2 x x +γ y 2 exp i 2π +ψ g(x, y) = exp −0.5 2σ 2 λ (6) where, λ is the wavelength of the sinusoidal factor, θ is the orientation, ψ is the phase offset, σ is the scale of the Gaussian envelope, γ is the spatial aspect ratio, x = x cos θ + y sin θ and y = −x sin θ + y cos θ . The Gabor filter response to the inverted green channel of the colored retinal image is obtained by a 2D convolution operator and is computed in the frequency domain. The detailed procedure can be seen in [46]. The maximum filter response over the angle θ , spanning [0, π ] in steps of π/18 is computed for each pixel in the image at different scales (σ = {3, 4, 5, 6}). The maximum response across the scales and orientation is taken as the value in the pixel feature vector. The Gabor filter response at scale 4 is illustrated in Fig. 5c.
123
Int J CARS
Image illumination and shade correction The retinal images in CHASE_DB1 contain the variation in background intensity due to non-uniform illumination. Therefore, in the pre-processing step, we aim to obtain a shade corrected image by removing the background lighting variations. For this purpose, the estimate of background I B E is obtained by applying a filtering operation with a large arithmetic mean kernel. The size of the filter kernel is not a critical parameter as long as it is large enough to ensure the blurred image contains no visible structures such as vessels. In this work, we use a 69×69 pixel size kernel. Then the difference between the morphologically opened image I O and the estimated background I B E is then computed for each pixel to obtain a background normalized image In . In (x, y) = I O (x, y) − IBE (x, y)
(7)
In this respect, there are shade-correction procedures reported in literature which are either based on subtracting the estimated background from original image [16,47] or on dividing the later by the former [48,49]. However, the experiments show [33] that the results of both of these methods are similar with no appreciable advantage of one on other. We have used the subtractive method in our work. Likewise, when the fluctuation in background intensity of retinal images is examined, there can be significant variation in intensities between images due to different illumination conditions in the acquisition process. Therefore, a shade corrected image is obtained by applying a global transformation with the purpose of reducing the intensity variation and contrast enhancement in relation to the green channel of RGB image. For this purpose the pixel intensities are modified according to the following global linear transformation function [33], ⎧ if In (x, y) < 0 ⎨ 0, if In (x, y) > 1 I H (x, y) = 1, ⎩ p(x, y), otherwise (8) p(x, y) = In (x, y) + 0.5 − gMAX_PIXELS where,I H (x, y) is the homogenized image, In (x, y) is the normalized image shown in (7), gMAX_PIXELS is the intensity value presenting the highest number of pixels in the normalized image In (x, y). The pixels with intensity value equal to gMAX_PIXELS belong to the background of the retinal image. This global transformation will set them to 0.5 and will standardize the intensity around this value of those background pixels with different illumination conditions, as illustrated in Fig. 6. Gradient orientation analysis The blood vessels are localized by analyzing the orientation of the gradient vector field. The unit gradient vectors of the
123
Fig. 6 a Retinal image, b normalized and shade corrected image
image are highly discontinuous along the bilaterally symmetrical regions, i.e., the linear structures which represent the blood vessels. Therefore, the blood vessels are localized by finding the discontinuities in the gradient orientation. The feature extraction depends on the orientation of the gradient vector field not its magnitude; therefore, it is robust against low contrast and non-uniform illumination. The gradient vectors for the image I (x, y) are approximated by the first-order derivative operators in the horizontal (kx) and vertical (ky) directions. gx (x, y) = I (x, y)∗ kx g y (x, y) = I (x, y)∗ ky
(9)
The gradient vectors gx(x, y) and gy(x, y) are normalized by dividing with their magnitude to compute the unit gradient vectors ux(x, y) and uy(x, y). u x (x, y) = gx (x, y)/ gx2 (x, y) + g 2y (x, y) (10) u y (x, y) = g y (x, y)/ gx2 (x, y) + g 2y (x, y) The unit vectors are assigned to zero if the gradient magnitude is too small (
ORIGINAL ARTICLE
Delineation of blood vessels in pediatric retinal images using decision trees-based ensemble classification Muhammad Moazam Fraz · Alicja R. Rudnicka · Christopher G. Owen · Sarah A. Barman
Received: 15 July 2013 / Accepted: 18 November 2013 © CARS 2013
Abstract Purpose Automatic segmentation of the retinal vasculature is a first step in computer-assisted diagnosis and treatment planning. The extraction of retinal vessels in pediatric retinal images is challenging because of comparatively wide arterioles with a light streak running longitudinally along the vessel’s center, the central vessel reflex. A new method for automatic segmentation was developed and tested. Method A supervised method for retinal vessel segmentation in the images of multi-ethnic school children was developed based on ensemble classifier of bootstrapped decision trees. A collection of dual Gaussian, second derivative of Gaussian and Gabor filters, along with the generalized multiscale line strength measure and morphological transformation is used to generate the feature vector. The feature vector encodes information to handle the normal vessels as well as the vessels with the central reflex. The methodology is evaluated on CHASE_DB1, a relatively new public retinal image database of multi-ethnic school children, which is a subset of retinal images from the Child Heart and Health Study in England (CHASE) dataset. Results The segmented retinal images from the CHASE_ DB1 database produced best case accuracy, sensitivity and specificity of 0.96, 0.74 and 0.98, respectively, and worst case measures of 0.94, 0.67 and 0.98, respectively. Conclusion A new retinal blood vessel segmentation algorithm was developed and tested with a shared database. The M. M. Fraz (B) · S. A. Barman School of Computing and Information Systems, Faculty of Science Engineering and Computing, Kingston University London, London, UK e-mail: [email protected] A. R. Rudnicka · C. G. Owen Division of Population Health Sciences and Education, St. George’s, University of London, London, UK
observed accuracy, speed, robustness and simplicity suggest that the algorithm may be a suitable tool for automated retinal image analysis in large population-based studies.
Keywords Medical imaging · Retinal image analysis · Ensemble classification · Image segmentation · Boot strapped decision tree
Introduction The human retina has the potential to reveal important information about retinal, ophthalmic, and even systemic diseases such as diabetes, hypertension, and arteriosclerosis [1]. Retinal imaging is increasingly used in large scale, populationbased studies for detection of glaucoma, diabetic retinopathy, age-related macular degeneration and cardiovascular disease [2]. A number of retinal blood vessel features such as arteriolar nicking and narrowing have been linked to systemic disease, and the morphological characteristics of retinal blood vessels themselves have been discussed in association with cardiovascular and coronary disease in adult life [3] and with retinopathy of prematurity in infancy [4]. The morphology of retinal vessels (particularly arterioles) has also been examined in relation to cardiovascular risk factors both in early and adult life [4]. Retinal vessels often have strong light reflexes along their centerline, which is more apparent on arteriolars than venules, and in younger compared with older participants, especially those with hypertension. The morphological characteristics of retinal images of premature infants and children are very different than that of the adult retina. Choroidal vessels are more visible alongside the retinal vessels in retinal images taken from premature infants [5]. A bright central reflex on the vessels and illumination artifacts
123
Int J CARS
Fig. 1 Vessel profiles: a–c Vessel profile of image from DRIVE dataset d–f Vessel profile of image from CHASE dataset: First column is the retinal image, second column is the magnified view of image part
marked in first column, and third column is the mesh plot of intensity profiles of second column
contribute to challenges in image processing when retinal images from school children are considered [4]. A variety of blood vessel segmentation algorithms is present in the literature [6]. However, most of them are evaluated on a limited range of retinal images of adult people, which includes twenty images each from the DRIVE [7] and the STARE [8] databases, which do not cater for the image-related characteristics such as, inter-image and intra-image variability in luminance, drift in contrast and uneven background gray level values, which are often found in the images acquired in larger population-based studies. The vessel profiles of a normal healthy adult retinal image and the image of a child retina from the CHASE [9] dataset are illustrated in Fig. 1, where the low contrast of blood vessels and the uneven background with augmented noise is clearly visible. Detecting blood vessels in retinal images of school children with the presence of uneven illumination, non-homogenized background, bright and dark lesions and central vessel reflexes is a challenging problem. This paper presents a supervised method for segmentation of blood vessels in retinal images of 9- and 10-yearold children of different ethnic origin using an ensemble classifier of bagged decision trees. The images named as CHASE_DB1 [9] are characterized with stark differences in background levels of retinal pigmentation (being more pigmented in South Asians and African Caribbean compared
with white Europeans). The classifier utilizes a feature vector based on the vessel-ness measure obtained by utilizing a group of the dual Gaussian matched filter, the Gabor filter and second-order derivative of Gaussian matched filter along with the multiscale line strength measure of blood vessels. The feature vector encodes information to successfully handle the blood vessels containing the central vessel reflex and without it. The classifier based on the boot strapped decision trees is a classic ensemble classifier which has been broadly applied in many application areas of image analysis, but has not been extensively utilized for retinal image analysis aimed at retinal vessel segmentation. The estimation of classification accuracy and prediction of feature importance during the training phase without using test data are the significant features of the boot strapped ensemble classifier. Moreover, the algorithm is computationally fast in training in addition to classification because it needs fewer samples for training as compared to other supervised algorithms available in the literature. The organization of the paper is as follows. In Sect. 2, the retinal vessel segmentation methodologies available in the literature are reviewed. The methodology, the feature vector and the classifier are explained in Sect. 3. In Sect. 4, the performance metrics, the accuracy and the robustness of the algorithm is presented and discussed. Finally, the paper is concluded in Sect. 5.
123
Int J CARS
Related work There is a substantial amount of work reported in the literature on automated segmentation of blood vessels in retinal images. A recent detailed review of these methods can be found in [6]. The vessel segmentation algorithms can be divided into two broad categories: rule-based methods and supervised methods. The rule-based methods can further be classified in to the techniques based on matched filtering, morphological processing, vessel tracking, multiscale analysis and model-based algorithms. The matched filtering methodology employs the piecewise linear approximation, the decrease in vessel diameter along vascular length and the Gaussian like intensity profile of retinal blood vessels and uses a kernel based on a Gaussian or its derivatives to enhance the vessel features in the retinal image. The matched filter was first proposed by Chaudhuri [10] and later adapted and extended by Hoover [11], Jiang [12] and Al-Rawi [13]. Mathematical morphology in combination with curvature evaluation [14] and matched filtering for centerline detection [15–17] is also exploited for retinal vessel segmentation. The algorithms based on vessel tracking [18,19] segment a vessel between two points using local information and work at the level of a single vessel rather than the entire vasculature. The multiscale approaches for vessel segmentation are based on scale-space analysis. The width of a vessel decreases gradually across its length as it travels radially outward from the optic disc. Therefore, the idea behind scale-space representation for vascular extraction is to separate out information related to the blood vessel having varying width. The multiscale second-order local structure of an image (Hessian) is examined and a vessel-ness measure is obtained on the basis of eigenvalue analysis of the Hessian [20]. A likelihood ratio test for vessel centerlines extraction that combines matched filter responses, confidence measures and vessel boundary measures was also presented [21]. The vessel cross-sectional intensity profiles can be approximated by a Gaussian or its higher-order derivatives, a mixture of Gaussians or hermite polynomial profile in cases where the central vessel reflex is present. Vermeer [22] modeled the vessel profile as a Laplacian kernel for segmentation. Mahadevan [23] presented a robust and modular framework for vessel detection in noisy images which is based on the estimation of the log likelihood of vessel parameters in a noisy environment. The framework is adaptable to incorporate a variety of vessel profile models including Gaussian, derivatives of Gaussian and dual Gaussian and various noise models like Gaussian noise and Poisson noise. A multiresolution hermite model was proposed for vascular segmentation by Li [24], which employed a two-dimensional hermite function intensity model in a quad-tree structure over a range of spatial resolutions. Lam et al. [25] proposed a vessel segmentation algorithm for pathological retinal images based on
the divergence of vector fields. The algorithm presented by Lam [26] was based on regularization- based multiconcavity modeling and was able to handle both normal and pathological retinas with bright and dark lesions simultaneously. Paulus et al. [27] proposed a methodology for quality assessment of retinal photographs. The cross-sectional intensity profile of retinal vessels was estimated by the dual Gaussian model initially used by Gao [28,29]. A number of authors have investigated the use of active contour models, or snakes [30] and geometric models based on level sets [31]. The supervised methods utilize ground truth data for the classification of vessels based on given features. These methods include the use of neural networks [32,33], principal component analysis [34], k-nearest neighbor classifiers [35], ridge-based primitives classification [36], Gaussian Mixture Model [37], line operators and support vector machine (SVM) classification [38], feature-based Ada-Boost classifier [39], SVM-based semi-supervised classification [40]. A decision tree-based ensemble classifier is presented in [41] for vessel segmentation which works well on normal and pathological retinal images. This work utilizes the ensemble classification approach aimed at segmenting the retinal vasculature from retinal images acquired from children attributed with normal vessels as well as the vessels with a central reflex.
The methodology The feature vector is designed to compute the quantifiable measure of the vessel-ness for each pixel in such a way that the ensemble classifier of bagged decision trees can successfully differentiate between the vessel (with and without the central vessel reflex) and the non-vessel region, an extension of work presented in [42]. A decision tree-based ensemble classifier is proposed in [41], where the main goal is to segment the blood vessels from adult retinal images containing pathologies. The color analysis of anatomical structures in retinal images reveals that the retinal pathologies appear either as dark lesions (e.g., hemorrhages) or bright lesions (soft and hard exudates). The utilized feature vector encodes information to successfully handle both normal and pathological retinas with bright and dark lesions simultaneously. In the current paper, we aim to solve the problem of retinal vessel segmentation from pediatric retinal images which are characterized with a bright central reflex in arterioles and wider venules; low contrast and uneven illumination. We have employed the feature vector which successfully differentiates the blood vessels containing the central vessel reflex and without it from the retinal background exhibiting poor contrast and uneven illumination. For this purpose the following changes are introduced in the feature vector.
123
Int J CARS
• A generalized multiscale line detector is used to compute the measure of the vessel-ness for each pixel in the retinal image. Besides highlighting the vessel structure, the generalized multiscale line detector reduces the false positive responses given to closely located vessels and at vessel crossovers, it deals with the central vessel reflex in an effective way. • The response of the Gaussian matched filter to the blood vessels with a central vessel reflex appears incorrectly as two blood vessels running parallel to each other. A twinGaussian model is used to give a correct response and imitate the vessels with a central light reflex which is apparent in the retinal images of children in the CHASE_DB1 database. • The retinal images in CHASE_DB1 contain a variation in background intensity due to non-uniform illumination. Therefore, we aim to obtain a shade corrected image by estimating and removing the background lighting variations using a simple linear transformation. The pixel intensity of the shade corrected image is included as one feature in the feature vector.
Moreover, we have also included a linear combination of line operator responses obtained at multiple scales, the morphological top-hat operation and gradient orientation analysis of the gradient vector field. The feature space is normalized to zero mean and unit standard deviation by applying a normal transformation, i.e., the filter response at each scale is subtracted by its mean and divided by its standard deviation. The inherent strength of ensemble classification for vessel segmentation is one in which multiple classifiers utilizing multiple cues/features for classification are strategically generated and combined to improve the classification or prediction performance of a model, or reduce the likelihood of a poor or unfortunate selection. This approach is intuitively used in our daily lives where we seek the guidance of multiple experts, weigh and combine their views in order to make a more informed and optimized decision. The boot strapped ensemble can estimate the classification accuracy and the feature importance index during the training phase without using the test data. The set of filters contain multiscale filters; therefore, the importance of the filter response at each scale is computed during classifier training, and the most significant features are selected to be included in the feature vector, thus reducing the dimensionality which in turn renders better computational performance. Moreover, the optimal number of decision trees used in ensemble creation can also be estimated during the training phase. The filters set and the line operator are applied to the inverted green plane of the colored RGB retinal image without pre-processing because the green channel of the RGB color fundus images gives the highest contrast between vessel and background. The green plane is
123
normalized and shade corrected in order to remove the variation in background illumination and shade correction before applying the morphological transformation and computing the Gradient Orientation Analysis (GOA) map. The feature vector Dual Gaussian and second derivative of Gaussian matched filter The small curvatures of blood vessels can be approximated by piecewise linear segments. The vessels appear in contrast relative to the background, have lower reflectance compared with the non-vessel retinal surface and almost never have ideal step edges. The intensity profile of the vessels can be approximated by a Gaussian curve and varies slightly from vessel to vessel. More specifically, the first-order derivative of Gaussian matched filter acts as the edge filter of the retinal vessel wall [6] and the second-order derivative of the Gaussian matched filter is employed for enhancement of piecewise linear structure of retinal vessel [43]. The Gaussian shaped matched filter and its extension have been used for vessel detection in retinal images [6]. The profile of the filter which is designed to match that of a blood vessel is defined as 1 x 2 exp − 2 ; for |x| f (x, y) = h 1 √ 2σ1 2π σ1 ≤ (t.σ1 ), |y| ≤ L/2
(1)
where, σ represents the scale of the filter; L” is the length of the neighborhood along the y-axis; h 1 is constant representing the height of the Gaussian curve, x = x cos θ + y sin θ, θ is the rotation angle of the filter kernel and t is a constant and is usually set to 3 because more than 99 % of the area under the Gaussian curve lies within the range [−3σ , 3σ ]. The response of the Gaussian matched filter to the blood vessels with a central vessel reflex looks like two blood vessels running parallel to each other. A twin-Gaussian model is used to imitate the vessels with a central light reflex which is apparent in the retinal images of children in the CHASE_DB1 database. For this purpose, a second small Gaussian curve is subtracted from the main Gaussian curve. The second Gaussian is oriented at the same angle and controlled by same parameters as the first one. The parameters of the second Gaussian control the height and width of central light reflex and the matched filter kernel is expressed as, 1 x 2 exp − 2 f (x, y) = h 1 √ 2σ1 2π σ1 1 x 2 exp − 2 − h2 √ (2) 2σ2 2π σ2
Int J CARS
Gaussian Matched Filter Kernel
Dual Gaussian Matched Filter Kernel
h1
h1
-X
0
+X
(a)
(b) Fig. 2 Dual Gaussian matched filter; a, b cross section and 2D profile of Gaussian matched filter, c, d cross section and 2D profile of dual Gaussian matched filter. (For the purpose of demonstration, the filter
-X
0
+X
(c)
(d) are constructed with following parameters; h 1 = 1, h 2 = 0.4∗ h 1 , σ1 = 2.5, σ2 = 0.8∗ σ1 , L = 27, t = 3 and = 0 degrees)
Fig. 3 Filter set responses; a retinal sub-image, b Dual Gaussian filter response highlighting the central reflex pixels, c second-order derivative of Gaussian matched filter response
In Fig. 2, the first column in the figure shows the cross section and two-dimensional plots of the Gaussian Matched filter kernel whereas the second column illustrates the cross section and two-dimensional profile of the dual Gaussian matched filter kernel. For the purpose of demonstration, the filter kernels shown in Fig. 2 are constructed with following parameters; h 1 = 1, h 2 = 0.4∗ h 1 , σ1 = 2.5, σ2 = 0.8∗ σ1 , L = 27, t = 3 and = 0 degrees. The second-order derivative of the Gaussian matched filter has also been exploited for vessel enhancement [43]. We have constructed a set of twenty four kernels in total,
consisting of two sets of twelve 2D dual Gaussian kernels and twelve second-order derivative of Gaussian kernels. The kernels in each set are oriented at twelve equiangular rotations θ1 · · · θ12 . The orientation angles of the kernels in each set of the filter group can be defined as {θ |0 ≤ θ ≤ π & θ mod(π/8) = 0}. The filter set is applied to the retinal image and the response is taken as the value for the highest scoring kernel at each pixel, which is taken as the pixel feature vector. Figure 3 illustrates a section of retinal image, the response of the dual Gaussian filter kernel and the secondorder Gaussian kernel. It is observed that the central vessel
123
Int J CARS
The retinal vasculature is composed of arteries and veins appearing as piecewise linear features, with variation in width and their tributaries visible within the retinal image. The concept of employing line operators for detection of linear structures in medical images is introduced in [44] which is modified and extended in [38,41,45] to incorporate the morphological attributes of retinal blood vessels. The average pixel intensity is evaluated along lines of a particular length passing through the pixel under consideration at 12 different orientations spaced by 15 degrees each. The line with the highest average pixel intensity is selected. The line strength of a pixel is calculated by computing the difference between the average gray values of a square sub-window centered at the target pixel with the average intensity of the selected line. A square sub-window of size w is centered at each pixel w . and average of pixel intensities is calculated, termed as Savg The average pixel intensity is measured along lines of a pixel length w passing through the center at 12 different orientations from 0 to 180 degrees at angular resolution of 15 degrees. The line with the highest average pixel intensity is w . The measure of the candidate line and is termed as Dmax vessel-ness of a pixel Vw is calculated by computing the difference between the average pixel intensity of a square subwindow centered at the target pixel with the average intensity of the candidate line as,
detector, the window size w is chosen in such a way that if the square sub-window is placed on the image and the center of this window corresponds to the center of vessel, then this window contains approximately the same number of vessel and background pixels. Therefore, it is selected to be twice the size of a typical vessel width in the image database. It has been experimentally found that W = 25 is a good choice for CHASE_DB1. The intensity of pixels in the middle of a vessel with a central reflex is usually lower than their surrounding pixels, instead of attaining the maximum value as in case of normal vessel. The line detector with longer lines can successfully recognize them as a vessel pixel because the candidate line includes only a small number of central reflex pixels; hence, the average intensity of pixels in the candidate line is not substantially affected and the central reflex pixels have a high vessel-ness measure. The vessel-ness image is shown in Fig. 4. It can be observed in Fig. 4c that the line detector with longer lines, i.e., a larger window size is effective in dealing with vessels with a central reflex but at the same time it produces false vessel responses at background pixels near the strongly contrasted vessels, i.e., those vessels which appear brighter than background. Also, it tends to merge the vessels in close proximity to each other and produces false responses in the area between vessels that are close to each other. These drawbacks can be avoided by using a generalized multi-scale line detector, which uses variable length of aligned lines in a fixed square sub-window. The vessel-ness measure using the generalized multi-scale line detector is can be expressed as;
w w − Savg Vw = Dmax
L w − Savg VwL = Dmax
reflex pixels are more highlighted with the dual Gaussian filter. Line strength features
(3)
The principal is that if the target pixel belongs to the vessel, then the vessel-ness measure will be large due to alignment of the candidate line along the direction of vessel. Contrary to that, the vessel-ness measure is low for a background pixel because the difference between the pixel intensity of the candidate line and square sub-window is small. In this basic line
(4)
where, 1 ≤ L ≤ w. The multi-scale line detector can be obtained by using different values of L. Besides highlighting the vessel structure, the shorter length line detector will not give false responses to closely located vessels and at vessel crossovers. Whereas, as described above, the longer length line detectors can deal with the central vessel reflex in an
Fig. 4 Line strength feature; a segment of retinal image; b, c Vessel-ness image with W = 3 and W = 21; d linear combination of multiscale line operator responses
123
Int J CARS
Fig. 5 Feature vector; a retinal image, b line strength feature, c Gabor filter response, d second derivative filter response, e gradient orientation analysis map, f morphological transformation
effective way. With the purpose of preserving the strength and exterminating the shortcoming of each individual line detector, the final vessel-ness measure is obtained by a linear combination of line responses computed by using line detectors at different scales. The combined vessel-ness measure can be defined as; 1 L Vw + I (5) Vcombined = nL + 1 L
where, n L is the number of scales used, VwL is the vessel-ness measure obtained using a line detector with line length L and square sub-window W , and I is the inverted green channel image. The combined vessel-ness image is shown in Figs. 4d and 5b, where it can be observed that a linear combination of responses from multi-scale line detectors works effectively in removal of background noise, dealing with the central vessel reflex while avoiding merging of close vessels and false responses near vessel crossovers. Multi-scale Gabor filter A Gabor filter is a linear filter and has been broadly used for multi-scale and multi-directional edge detection. The Gabor
filter can be fine-tuned to particular frequencies, scales and directions and therefore acts as low level feature extractor and background noise suppresser. The impulse response of a Gabor filter kernel is defined by the product of a Gaussian kernel and a complex sinusoid. It can be expressed as
2 x x +γ y 2 exp i 2π +ψ g(x, y) = exp −0.5 2σ 2 λ (6) where, λ is the wavelength of the sinusoidal factor, θ is the orientation, ψ is the phase offset, σ is the scale of the Gaussian envelope, γ is the spatial aspect ratio, x = x cos θ + y sin θ and y = −x sin θ + y cos θ . The Gabor filter response to the inverted green channel of the colored retinal image is obtained by a 2D convolution operator and is computed in the frequency domain. The detailed procedure can be seen in [46]. The maximum filter response over the angle θ , spanning [0, π ] in steps of π/18 is computed for each pixel in the image at different scales (σ = {3, 4, 5, 6}). The maximum response across the scales and orientation is taken as the value in the pixel feature vector. The Gabor filter response at scale 4 is illustrated in Fig. 5c.
123
Int J CARS
Image illumination and shade correction The retinal images in CHASE_DB1 contain the variation in background intensity due to non-uniform illumination. Therefore, in the pre-processing step, we aim to obtain a shade corrected image by removing the background lighting variations. For this purpose, the estimate of background I B E is obtained by applying a filtering operation with a large arithmetic mean kernel. The size of the filter kernel is not a critical parameter as long as it is large enough to ensure the blurred image contains no visible structures such as vessels. In this work, we use a 69×69 pixel size kernel. Then the difference between the morphologically opened image I O and the estimated background I B E is then computed for each pixel to obtain a background normalized image In . In (x, y) = I O (x, y) − IBE (x, y)
(7)
In this respect, there are shade-correction procedures reported in literature which are either based on subtracting the estimated background from original image [16,47] or on dividing the later by the former [48,49]. However, the experiments show [33] that the results of both of these methods are similar with no appreciable advantage of one on other. We have used the subtractive method in our work. Likewise, when the fluctuation in background intensity of retinal images is examined, there can be significant variation in intensities between images due to different illumination conditions in the acquisition process. Therefore, a shade corrected image is obtained by applying a global transformation with the purpose of reducing the intensity variation and contrast enhancement in relation to the green channel of RGB image. For this purpose the pixel intensities are modified according to the following global linear transformation function [33], ⎧ if In (x, y) < 0 ⎨ 0, if In (x, y) > 1 I H (x, y) = 1, ⎩ p(x, y), otherwise (8) p(x, y) = In (x, y) + 0.5 − gMAX_PIXELS where,I H (x, y) is the homogenized image, In (x, y) is the normalized image shown in (7), gMAX_PIXELS is the intensity value presenting the highest number of pixels in the normalized image In (x, y). The pixels with intensity value equal to gMAX_PIXELS belong to the background of the retinal image. This global transformation will set them to 0.5 and will standardize the intensity around this value of those background pixels with different illumination conditions, as illustrated in Fig. 6. Gradient orientation analysis The blood vessels are localized by analyzing the orientation of the gradient vector field. The unit gradient vectors of the
123
Fig. 6 a Retinal image, b normalized and shade corrected image
image are highly discontinuous along the bilaterally symmetrical regions, i.e., the linear structures which represent the blood vessels. Therefore, the blood vessels are localized by finding the discontinuities in the gradient orientation. The feature extraction depends on the orientation of the gradient vector field not its magnitude; therefore, it is robust against low contrast and non-uniform illumination. The gradient vectors for the image I (x, y) are approximated by the first-order derivative operators in the horizontal (kx) and vertical (ky) directions. gx (x, y) = I (x, y)∗ kx g y (x, y) = I (x, y)∗ ky
(9)
The gradient vectors gx(x, y) and gy(x, y) are normalized by dividing with their magnitude to compute the unit gradient vectors ux(x, y) and uy(x, y). u x (x, y) = gx (x, y)/ gx2 (x, y) + g 2y (x, y) (10) u y (x, y) = g y (x, y)/ gx2 (x, y) + g 2y (x, y) The unit vectors are assigned to zero if the gradient magnitude is too small (
