Computerized Medical Imaging and Graphics 37 (2013) 337–345
Contents lists available at ScienceDirect
Computerized Medical Imaging and Graphics journal homepage: www.elsevier.com/locate/compmedimag
Reliable monitoring system for arteriovenous ratio computation S.G. Vázquez a,∗ , N. Barreira a , M.G. Penedo a , M. Rodríguez-Blanco b a b
Varpa Group, Department of Computer Science, University of A Coru˜ na, Spain Service of Ophthalmology, Complejo Hospitalario Universitario, Santiago de Compostela, Spain
a r t i c l e
i n f o
Article history: Received 27 November 2012 Received in revised form 10 September 2013 Accepted 7 October 2013 Keywords: Retinal analysis Arteriovenous ratio AVR computation Retinal vessel classification AVR monitoring system
a b s t r a c t The degree of narrowing or widening in retinal vessels related to several cardiovascular diseases such as hypertension or diabetes may be measured by the arteriovenous ratio (AVR), that is, the relation between the artery and vein retinal vessel widths. Nevertheless, its lack of reproducibility, due mainly to a laborious manual calculation and the dependence of the vessels selected for its estimation, hinders its use in daily medical practice. This variation makes difficult to monitor the patient’s condition over time. This paper describes a reliable AVR monitoring system which computes automatically the AVR from several images of the same patient acquired at different times using the same vessels measured at the same points. The system has been evaluated in a large data set of 158 pairs of images and good correlation results between medical experts and the system have been achieved. © 2013 Elsevier Ltd. All rights reserved.
1. Introduction Retinal blood vessels constitute the only vascular network in the human body which can be directly visualized non-invasively. Thus, the increasing development of image analysis techniques in the retinal field provides researchers and physicians the means to study the relationship between retinal micro-circulation and macro-circulation. And therefore, they constitute the key for an advanced diagnosis of certain diseases which cause early changes in retinal vascular network, such as tortuosity or alterations in vessel widths. Thus, arteriolar and venular changes seem to reflect different pathologies. Arteriolar narrowing is considered an initial damage in hypertension, whereas venular widening is associated with hyperglycemia and obesity. Moreover, they are both associated with stroke and coronary heart disease [1]. Hence, the retinal vascular width is a fundamental parameter in the study and early diagnosis of these pathologies. In this manner, the arteriovenous ratio (AVR), that is, the relation between arteriolar and venular vessel widths, is a popular dimensionless measure to assess the patient’s condition. The AVR is computed as the ratio between artery and vein widths measured in several circumferences centered at the optic disc. In the literature, several methods have been proposed to estimate this measure. In [2,3], the AVR is computed as the quotient
∗ Corresponding author. Tel.: +34 981167000. E-mail addresses: [email protected] (S.G. Vázquez), [email protected] (N. Barreira), [email protected] (M.G. Penedo), [email protected] (M. Rodríguez-Blanco). 0895-6111/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compmedimag.2013.10.001
between the averages of several arteriolar and venular widths. Benavent et al. consider that the AVR is difficult to interpret and therefore they proposed using the average of artery and vein widths separately [4]. Parr and Spears [5,6] and subsequently Hubbard et al. [7] derived formulas widely used in popular studies and other methods [8]. In these methods, the AVR is computed as the quotient of two variables, the central retinal artery equivalent (CRAE) and the central retinal vein equivalent (CRVE), which represent the relation among a vessel trunk and its two branches. The formulas to estimate the equivalents have been derived theoretically and empirically, observing a fixed data set, and some reformulations have been introduced. In fact, Knudtson et al. [9] realized that the Hubbard’s formulas were dependent on the number of selected vessels and they propounded a reformulation using only the six main arteries and veins. Recently, Patton et al. [10] proposed a revised CRAE formula, where the branching coefficient, BC, is not constant as in the Knudtson’s formula (BC = 1.28) but a linear function dependent on the asymmetry index, AI, BC = 0.78 + 0.63 * AI, where the asymmetry index is the quotient between the width of the two branches. The previous methods still require user assistance in some phases of the methodology, such as the distinction between arteries and veins, the optic nerve location or the selection of vessels for the AVR computation. Recently, some automatic approaches have been presented [11–14]. Despite the previous methods provide an AVR measure valid for research purposes, its applicability to daily clinical practice has not been established. In fact, the AVR is affected, in general, by imprecision and subjectivity caused by several sources of variation, such as the camera influence or the manner to compute the
338
S.G. Vázquez et al. / Computerized Medical Imaging and Graphics 37 (2013) 337–345
value [1]. In general, as Knudtson realized, the AVR value depends on the selected set of vessels and the area where they are measured. Vazquez et al. [14] show this influence in case of the computation of the AVR as an average width ratio. But the same issue could happen using the Knudtson’s formulas when, for example, it is not possible to detect six veins and arteries in the image. In addition, in several studies, the AVR values estimated by an expert and an automatic system are similar in average but the correlation between them is low, even, when the correlation was calculated between the AVR values achieved by medical experts using semiautomatic tools [13,14]. Thus, for all of the above reasons, reliable AVR computations should involve the measurement of the same set of vessels at the same points for each patient’s eye. In this work, we propose a two-stage methodology to compute the patient’s AVR over time in a reliable and repeatable manner. In the first stage, the algorithm calculates the AVR of a sample image taken at time t0 . For that, it locates the optic nerve and selects several radii of interest where the vessels are detected and measured. After that, the vessel segments are classified into veins and arteries and the AVR value is computed from a subset of the detected vessel segments [14]. The second stage allows to compute automatically the AVR in new sample images of the same patient’s eye acquired at times ti by registration. The image taken at time t0 is considered as the reference image and the AVR measurements in subsequent images are estimated in set of vessel segments selected for the reference image. The proposed methodology has been implemented as an AVR monitoring service on the SIRIUS system [15]. This system is a webbased application for retinal image analysis and patient’s checkup management. In this manner, we can assess the patient’s condition over time. The SIRIUS system also provides an automatic and a semiautomatic services to compute the AVR. At the semiautomatic mode, the system detects the vascular tree, measures the vessel widths and estimates the AVR. However, it still requires the user input to locate the optic nerve, select the vessels suitable for the calculus and classify them into arteries and veins. The paper is structured as follows. Section 2 explains the AVR monitoring method. In Section 3 the data set used is presented. Section 4 describes the performed experiments and shows the achieved results. Finally, in Section 5, we provide some concluding remarks. 2. Method The proposed method to estimate the AVR is based on a vessel registration approach which allow to measure the vessel widths at the same points from different sample images taken at different dates. This method consists of two main stages. In the first stage, a reference AVR is estimated using the application SIRIUS in a semior automatic way. This result is used to achieve subsequent AVR results by a vessel registration approach. The sample image taken at time t0 where the reference AVR is computed constitutes the reference image. In the second stage, we compute the AVR automatically in a new sample image from the same patient taken at time ti . First, this new image is registered to the reference image. Then, we locate in the new image the vessel segments in the reference image and we measure their widths. Finally, we apply the selected formula to get the AVR value. Fig. 1 shows a schema of the monitoring process. Fig. 2 includes two screenshots which show the integration of the AVR monitoring system in the SIRIUS web application. 2.1. Reference AVR measurement The reference AVR is measured by a module of the SIRIUS system [15]. The system has two modes of operation to estimate the
AVR. The first mode is semiautomatic and requires the user input to locate the optic disc, select the set of suitable vessels for the calculus and classify them into arteries and veins. The second mode is fully automatic, however, the user can supervise the result, discarding vessels, considering new ones or modifying the classification given by the system. In this section, we briefly describe the stages of the proposed automatic methodology. This comprises several stages. First, the optic nerve is located and then the vessels are segmented and measured in several circumferences concentric to the optic disc. After that, the vessel segments obtained are classified as arteries and veins. Finally, the system selects an appropriate set of vessel segments to estimate the AVR. 2.1.1. Optic disc location The optic disc location is based on the Hough transform [16]. First, we delimit the region of interest by a blob detection method, the Difference of Gaussians operator. For that, the original image is convolved with two Gaussian filters at different scales and the difference between these convolutions is obtained. The maxima of the DoG output corresponds to a rough location of the optic disc center. Then, we crop the image around this point using domain information and we apply the Hough transform. In order to restrict the set of points that can vote in the Hough transform, we detect the edges with the Sobel operator, whereas a crease detector algorithm is applied to obtain the vessel centerlines. If we consider the image as a topographic relief, where the pixel intensities represent the surface elevation, a crease is defined as the set of continuous image points that links the highest or lowest elevation levels. The crease detector used was the Multilocal Level Set Extrinsic Curvature enhanced by the Structure Tensor (MLSEC-ST) operator proposed in [17]. Thus, combining the Sobel and MLSEC-ST outputs, we obtain the optic disc perimeter points, these are, points that belong to edges but which are not near a vessel centerline. Once we have some optic disc perimeter points, we extract the optic nerve using the Hough transform for circle detection. 2.1.2. Measurement of vessel widths Once the optic disc is located, the vessels widths are measured in several circumferences centered on the optic disc. The number and radii of the circumferences can be changed by users in the semiautomatic mode. In the automatic mode, the circumferences are equally spaced a value of 0.5DD, being DD the optic disc diameter, and their radii range from 1DD to the border of the field-of-view (FOV). The proposed method to measure the vessel widths is based on snakes and it was originally proposed in [18]. A snake is a contour defined within the image which evolves guided by internal and external forces to fit some features. At first, the vessel centerlines are detected by the MLSEC-ST operator [17]. Thus, for a circumference of radius r, each snake consist of two chains of nodes initialized at the intersection points between each vessel centerline and the circumferences of radii r − n and r + n, with n = 10 pixels. Both chains of nodes advance along the normal to the vessel centerline in opposite directions toward the vessel walls guided by inner and outer energies. There are two types of nodes: common nodes which fit the edges and corner nodes which are adjusted to the corners of the vessel segment to be measured. In this model, the internal energy can be written as:
2 ∂s (s) ∂ s (s) Eint ((s)) = ˛(s) + ˇ(s) ∂s ∂ s2
(1)
where s (s) = [x(s), y(s)] defines the snake with node coordinates x(s), y(s), whereas the first and second term represent the first and second order derivatives, respectively. The parameters ˛(s) and ˇ(s)
S.G. Vázquez et al. / Computerized Medical Imaging and Graphics 37 (2013) 337–345
339
Fig. 1. Schema of the AVR monitoring process.
controls the snake shape. These values are constant for common nodes, in particular, ˛(s) = 0.25 and ˇ(s) = 0.01. The internal energy is null for corner nodes (˛(s) = ˇ(s) = 0) which implies the nodes are only attracted by the external forces. The external energy which attracts the snake to the target edge is defined by the next equation: Eext ((s)) = + ıEdist ((s)) + Egrad ((s)) + ωEstat ((s))
(2)
where is the dilation pressure, a vectorial magnitude which adjusts for each snake node its direction and sense of advance, Edist (i ) is the euclidean distance between a node and the vessel edges computed by the Canny operator from the input gray level image. These two terms are responsible for reaching the edges. On the contrary, the latter terms are stopping forces. Thus, the gradient energy, Egrad ((s)) brakes a node if its gradient is negative, that is, contrary to the advance direction. Finally, the stationary energy is based on the idea that probably a node should stop (keep on moving) if their adjacent nodes have stopped (kept on moving). This energy prevents movements of a node when its neighboring nodes have stopped or that a node stops in a local minimum when its neighboring nodes keep on moving. After the energy minimization, the snake contour should have a parallelogram shape and the nodes should be positioned in the vessel edges. However, image noise and discontinuities in the edge image may cause some nodes are placed outside or inside the edges. Hence an adjustment of snake contour is needed to correct the wrong nodes. Thus, a simple linear regression model is applied to the nodes at each snake side, achieving two regression lines. Then, the nodes are reorganized to fit to the corresponding regression line. If the correction is not possible, the snake is ruled out since it is inappropriate to be measured. In addition, if the final snake configuration has less than two nodes at each side corresponding to a vessel edge because many nodes evolved to the same position, the snake is also ruled out. Finally, the vessel width is estimated taken into account seven profiles in the snake parallelogram perpendicular to the centerline and equally spaced. At each profile, two distance measures are computed: the euclidean distance between the profile end points,
ω ˆ euc , and the distance obtained from a parabolic regression model, ω ˆ p . This last measure is used to correct the edge shift caused by a Gaussian smoothing applied to remove noise before applying the Canny edge detector. Thus, the profile gray levels are fitted to a parabola. To this end, only five points in the profile are considered: the end points, the point over the vessel centerline and the two points on either side of the latter with maximum gray level. Hence, the width of each profile is computed as the average of both measures (ω ˆ = (ω ˆ euc + ω ˆ p )/2). The final vessel width is estimated by the average of the three central values of the seven profile widths. In the automatic mode, before the vessel width measurement, the bifurcations and crossovers are detected [19] and they are ruled out. This step is necessary because both the measurement and classification of these points may introduce errors in the final AVR. 2.1.3. Vessel classification into artery and vein The classification of retinal vessels into arteries and veins is not straightforward. The biggest difference between arteries and veins is the color since veins are darker than arteries. However, the significance of this difference depends on the vessel size, being imperceptible for small vessels. In addition, the inter and intraimage lightness variability make the task more difficult. In [20], we have proposed an unsupervised method to classify the vessels which is robust to lightness and color variations. At each circumference, a coordinate axis centered on the optic disc divides the retina image into four quadrants. At each quadrant the vessel segments detected in the previous step are classified using the k-means clustering algorithm. Then, the coordinate axis is rotated through an angle of 20◦ and a new classification procedure is performed at each new quadrant. These two steps are repeated until it is achieved a 180◦ rotation. Then, the final class of a vessel segment is obtained, combining the results achieved in all the quadrants where the segment was classified. In order to avoid the influence of color variation, the input image is normalized using the multi-retinex technique [21]. Finally, we apply a procedure based on the minimal paths approach [22] to join vessel segments in consecutive circumferences which belong to the same vessel. Thus, the local results
340
S.G. Vázquez et al. / Computerized Medical Imaging and Graphics 37 (2013) 337–345
Fig. 2. Integration of the AVR monitoring system in the SIRIUS web application. (a) Interface for visualization and edition the AVR in the new sample image (first tab) and the reference image (second tab). (b) Interface which shows at right the list of the AVR measurements computed in new sample images using the reference AVR displayed at left side.
obtained in all the connected circumferences are combined to ensure the classification.
2.1.4. AVR computation Once the vessel segments have been classified and measured, we estimate the AVR as the ratio between the average artery and vein vessel widths [3,2]. Nevertheless, in the automatic mode, we make a vessel selection previously. In addition to the crossovers and bifurcations, we eliminate from the calculus the vessel segments whose width is over- or sub-estimated in relation with the rest of their connected vessel segments. According to experts, we do not take into account thin vessels, as well. Thus, we remove the arteries and veins whose widths are smaller than the 10-percentile [14].
Once the reference AVR has been estimated and saved for future computations, we also store some useful data for the vessel registration procedure, in particular, the optic disc center, the class of each vessel segment and the intersection points between each analysis circumference and the vessel centerline where the snakes have been initialized.
2.2. AVR measurement in a new sample image Given a new sample image taken at time ti , the method computes its AVR by registration using the stored reference AVR. For that, first, an image registration process between the two images is performed. Then the correspondence of the reference vessel
S.G. Vázquez et al. / Computerized Medical Imaging and Graphics 37 (2013) 337–345
341
Table 1 Image registration accuracy in a dataset of 20 retinal images.
Mean Standard deviation
a − r
MSE
RA(%)
0.033 0.050
0.097 0.070
99.100 1.550
and the number of seeds considered at each level were 6, 3 and 1, respectively. In the remaining levels, a Downhill Simplex algorithm is applied to optimize the correlation. At those levels, the search seeds are the optimized transformations of the previous level and the search process finishes when the difference between the maximum and the minimum values found in a neighborhood is lower than a preestablished threshold. Once the registration process has been finished and the two images are aligned, the registered crease images are used to obtain a similarity measure between them in order to determine the quality of the register process. This measure is the normalized cross-correlation coefficient, defined as:
Fig. 3. Schema of the multiple resolution registration process.
segments in the new sample image is obtained in order to measure their widths in the new image and compute the AVR in the same points. 2.2.1. Image registration The image registration algorithm used is a feature-based affine registration method whose landmark is the vessel tree [23]. Following the idea that the vessels can be thought as creases (ridges and valleys) when images are seen as landscapes, the registration process is based on the alignment of crease images. Thus, the first stage of the algorithm is to obtain the crease images of the reference and the new sample images using the MLSEC-ST operator [17]. After obtaining the crease images, an iterative optimization process to align both images is performed. That is, the reference image is fixed and the other one is transformed until a global maximum is accomplished. A suitable function to measure the alignment quality is the correlation function: Corr =
f (x)g((x))
(3)
x∈f
where f and g are the crease images and represent the transformation whose five parameters (x and y translation, tx , ty , rotation angle in clockwise direction and x and y scale, sx , sy ) we want to test. The search space is defined by the function Corr and the five parameters of the transformation. The function is non monotonic, that is, it has many local maxima and it is expensive to compute since it involves image transformations. Thus, to simplify the optimization process, the multiple resolution approach proposed by Elsen et al. [24] is carried out. Moreover, with the same aim in mind, the transformation is not applied to all pixels in the image, only the main creases with pixel values higher than a fixed threshold are transformed. The multiple resolution approach is handled by two pyramids (see Fig. 3) where the two crease images to register are at the bottom and each pyramid level is a half-resolution version of the image in the previous level, until images have a size of 64 pixels in each dimension. Then, the search starts at the top of the pyramid where an exhaustive search is accomplished and the correlation is computed in the Fourier domain. Thus, the search seeds are the image transformations which maximize Corr computed in the frequency domain. At this level, the transformations performed involve only rotations of an angle of 5◦ during a certain number of iterations. For a pyramid of three levels, the number of iterations was 10
=
x,y
x,y
[f (x, y) − f ] [g(x, y) − g ]
[f (x, y) − f ]
2
x,y
[g(x, y) − g ]
2
(4)
where g is the mean of the registered image and f is the mean of the reference image. The registration process is valid if is higher than a threshold, otherwise, the AVR cannot be computed using monitoring. In order to evaluate the registration method, we use a dataset of 20 retinal images taken with a Canon CR6-45NM non-mydriatic retinal camera with a 768 × 584 pixel resolution. For each image, 50 random transformations were applied with maximum translation and rotation values of ±100 pixels and ±5◦ , respectively, and without scaling. These transformations represent the maximum transformation values of the majority of the analyzed images. Table 1 shows the mean and standard deviation values of the absolute difference between the transformation applied and the transformation recovered by the method, the mean square error (MSE) of the recovered transformations and the registration accuracy (RA). The absolute difference was computed as a − r = |(tx − tx ) + (ty − ty ) + ( − ) + (sx − sx ) + (sy − sy )|, where a and r represent the applied and the recovered transformations being tx , ty , , sx , sy and tx , ty , , sx , sy their translation, rotation and scaling parameters, respectively. The registration accuracy represents the percentage of the transformations where the applied and the recovered transformations are equal with ±1 pixel of difference. To obtain this measure, four points are selected manually and their positions in the transformed and the dynamic images are compared. Thus, the average of the euclidean distances between the positions of the four points is calculated. If the absolute value of the average is less or equal than 1 pixel, the applied and the recovered transformation are aligned. The small MSE value and the high RA indicate the good performance of the image registration method. 2.2.2. Vessel registration Once a suitable transformation matrix between the reference image and new sample image has been obtained in the previous step and the transformation quality has been proved, we can measure the vessel widths in the new image at the same points considered in the reference image. First, the stored reference optic disc and intersection points are transformed into the new sample image (Fig. 4). For each transformed intersection point, we measure its distance, d, to the transformed optic disc center and we measure the width of the corresponding vessel segment by a snake. Once the snake evolution finishes, the topological adjustment is applied as well, and the vessel width is measured. The classification of the
342
S.G. Vázquez et al. / Computerized Medical Imaging and Graphics 37 (2013) 337–345
Fig. 4. Schema of the vessel registration procedure. In the new image (right), the vessel widths are measured at the intersection points registered from the reference image (left). Note that the same set of vessel segments is used to estimate the AVR.
vessels segments in the new image is obtained from the corresponding vessel segment in the reference image. Finally, the AVR is computed as the ratio between the average artery and vein vessel widths. Note that, we use the middle point of the snake initialization for the vessel registration and not the adjusted snake points. This is due to the fact that the final snake could be quite away from the initial position because of its deformation. If the snakes evolve from the same starting points in both images, there is a higher probability of similar deformations.
less than 2.2e−16 and a Pearson’s correlation coefficient of 0.61. Moreover, a p-value of 0.20 obtained in the second test reveals the AVR values are equal in average. However, the correlation coefficient and the p-value for the mean difference test increase to 0.92 and 0.89, respectively in a data set of 90 images, ruling out those images with the highest differences between the AVR values (see Table 2). The dispersion graphs of Fig. 5 show the correlation results. If we analyze the 25 images which present the AVR difference between expert and system highest (AVR6MMonitoring − AVR6MExpert ≥0.1) we derive the following assumptions:
3. Materials Eighty-nine hypertensive patients were studied from the Conxo Hospital in Santiago de Compostela (Spain). Two images, one per eye, have been acquired with a Cannon CR6-45NM non-mydriatic retinal camera for almost all patients. The images were stored with a resolution of 768 × 576 pixels in JPEG format. Each patient, after a basal assessment and an image acquisition, started an anti-hypertensive treatment. After 6 months, a new image was acquired. An ophthalmologist with many years of experience in the study of hypertensive retinopathy has analyzed the 158 pairs of images taken before and after treatment. The arteriovenous ratio in this set of images was computed by the expert using the SIRIUS system in a semiautomatic mode. In order to evaluate the monitoring system proposed in this work, each image taken at month 0 is considered the reference image and its AVR value graded by the expert represents the AVR reference. From the AVR reference, we compute the AVR in the image taken at month 6 and we compare the result obtained with the AVR graded by the expert in this image. 4. Results The proposed methodology has been tested in the data set composed of 158 pairs of images. We compute the AVR in the images taken at month 6 using the corresponding image captured at month 0 as reference image and its AVR value graded by the expert as AVR reference. A statistical analysis has been performed to test if the monitoring system behaves like an human expert, that is, whether despite using different vessels for the estimation, the AVR values obtained by the expert and the monitoring system are similar. We denote by AVR6MMonitoring the AVR obtained by the monitoring system and with AVR6MExpert the value measured by the medical expert in the images taken at month 6. A two-sided hypothesis test on Pearson’s correlation coefficient (H0 : AVR6MMonitoring ,AVR6MExpert = 0), as well as, a two-sided Welch’s t-test (H0 : AVR6MMonitoring − AVR6MExpert = 0) have been conducted. The first test provides evidence of correlation by means of a p-value
• The number of vessels taken into account to compute the AVR differs greatly from the image graded by the expert to the image measured by the monitoring system. This occurs mainly because the expert has selected a different number of vessels at months 0 and 6 as Fig. 6 shows or because the system does not detect some vessels at month 6 due to the image quality. As experts use a small number of vessels, these differences greatly affect the AVR value since it is computed as the ratio between the averages. • In other cases, it seems that the AVR value is different because the expert has selected in only one of the images a very characteristic vessel (wide or narrow) with a huge influence in the AVR (Fig. 6). • Moreover, an incorrect snake deformation is also possible, resulting in an under- or overestimation of the vessel width in one image whereas the width was properly measured in the other image. 4.1. Improvement of monitoring system First, we try to improve the system by preventing the elimination of vessels. The registration process works in the approved manner, so the vessel detection and measurement stages are the responsible of the elimination. Hence, we modify the measurement method, so each time a snake is invalidated due to its small number of nodes, it is thrown again, using as seed point a pixel located in an 8-neighborhood of the initial seed. The snake is definitely invalidated if any of the 8 possible seeds produces an enclosed snake. After introducing this change, the Pearson’s correlation coefficient in the 158 images increases slightly to 0.6329 while the p-value of the mean difference test remains the same, 0.209. Most of the missing vessels in the new sample image are lost because they are not detected. Their centerlines are not detected because the image quality. Thus, we modify the parameters of the centerline operator, which let to detect a greater number of vessels. These parameters are optimized to avoid false positives, but, since the vessel positions are known in the sample image, we can be less restrictive. Despite the change, the correlation value (0.6325) did not improve. This could happen because these parameters have
S.G. Vázquez et al. / Computerized Medical Imaging and Graphics 37 (2013) 337–345
343
Table 2 Test hypothesis results between the AVR values obtained by the human expert and the system in the images taken after 6 months of treatment. Pearson’s test p-value
Welch’s test p-value
0.616
Contents lists available at ScienceDirect
Computerized Medical Imaging and Graphics journal homepage: www.elsevier.com/locate/compmedimag
Reliable monitoring system for arteriovenous ratio computation S.G. Vázquez a,∗ , N. Barreira a , M.G. Penedo a , M. Rodríguez-Blanco b a b
Varpa Group, Department of Computer Science, University of A Coru˜ na, Spain Service of Ophthalmology, Complejo Hospitalario Universitario, Santiago de Compostela, Spain
a r t i c l e
i n f o
Article history: Received 27 November 2012 Received in revised form 10 September 2013 Accepted 7 October 2013 Keywords: Retinal analysis Arteriovenous ratio AVR computation Retinal vessel classification AVR monitoring system
a b s t r a c t The degree of narrowing or widening in retinal vessels related to several cardiovascular diseases such as hypertension or diabetes may be measured by the arteriovenous ratio (AVR), that is, the relation between the artery and vein retinal vessel widths. Nevertheless, its lack of reproducibility, due mainly to a laborious manual calculation and the dependence of the vessels selected for its estimation, hinders its use in daily medical practice. This variation makes difficult to monitor the patient’s condition over time. This paper describes a reliable AVR monitoring system which computes automatically the AVR from several images of the same patient acquired at different times using the same vessels measured at the same points. The system has been evaluated in a large data set of 158 pairs of images and good correlation results between medical experts and the system have been achieved. © 2013 Elsevier Ltd. All rights reserved.
1. Introduction Retinal blood vessels constitute the only vascular network in the human body which can be directly visualized non-invasively. Thus, the increasing development of image analysis techniques in the retinal field provides researchers and physicians the means to study the relationship between retinal micro-circulation and macro-circulation. And therefore, they constitute the key for an advanced diagnosis of certain diseases which cause early changes in retinal vascular network, such as tortuosity or alterations in vessel widths. Thus, arteriolar and venular changes seem to reflect different pathologies. Arteriolar narrowing is considered an initial damage in hypertension, whereas venular widening is associated with hyperglycemia and obesity. Moreover, they are both associated with stroke and coronary heart disease [1]. Hence, the retinal vascular width is a fundamental parameter in the study and early diagnosis of these pathologies. In this manner, the arteriovenous ratio (AVR), that is, the relation between arteriolar and venular vessel widths, is a popular dimensionless measure to assess the patient’s condition. The AVR is computed as the ratio between artery and vein widths measured in several circumferences centered at the optic disc. In the literature, several methods have been proposed to estimate this measure. In [2,3], the AVR is computed as the quotient
∗ Corresponding author. Tel.: +34 981167000. E-mail addresses: [email protected] (S.G. Vázquez), [email protected] (N. Barreira), [email protected] (M.G. Penedo), [email protected] (M. Rodríguez-Blanco). 0895-6111/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compmedimag.2013.10.001
between the averages of several arteriolar and venular widths. Benavent et al. consider that the AVR is difficult to interpret and therefore they proposed using the average of artery and vein widths separately [4]. Parr and Spears [5,6] and subsequently Hubbard et al. [7] derived formulas widely used in popular studies and other methods [8]. In these methods, the AVR is computed as the quotient of two variables, the central retinal artery equivalent (CRAE) and the central retinal vein equivalent (CRVE), which represent the relation among a vessel trunk and its two branches. The formulas to estimate the equivalents have been derived theoretically and empirically, observing a fixed data set, and some reformulations have been introduced. In fact, Knudtson et al. [9] realized that the Hubbard’s formulas were dependent on the number of selected vessels and they propounded a reformulation using only the six main arteries and veins. Recently, Patton et al. [10] proposed a revised CRAE formula, where the branching coefficient, BC, is not constant as in the Knudtson’s formula (BC = 1.28) but a linear function dependent on the asymmetry index, AI, BC = 0.78 + 0.63 * AI, where the asymmetry index is the quotient between the width of the two branches. The previous methods still require user assistance in some phases of the methodology, such as the distinction between arteries and veins, the optic nerve location or the selection of vessels for the AVR computation. Recently, some automatic approaches have been presented [11–14]. Despite the previous methods provide an AVR measure valid for research purposes, its applicability to daily clinical practice has not been established. In fact, the AVR is affected, in general, by imprecision and subjectivity caused by several sources of variation, such as the camera influence or the manner to compute the
338
S.G. Vázquez et al. / Computerized Medical Imaging and Graphics 37 (2013) 337–345
value [1]. In general, as Knudtson realized, the AVR value depends on the selected set of vessels and the area where they are measured. Vazquez et al. [14] show this influence in case of the computation of the AVR as an average width ratio. But the same issue could happen using the Knudtson’s formulas when, for example, it is not possible to detect six veins and arteries in the image. In addition, in several studies, the AVR values estimated by an expert and an automatic system are similar in average but the correlation between them is low, even, when the correlation was calculated between the AVR values achieved by medical experts using semiautomatic tools [13,14]. Thus, for all of the above reasons, reliable AVR computations should involve the measurement of the same set of vessels at the same points for each patient’s eye. In this work, we propose a two-stage methodology to compute the patient’s AVR over time in a reliable and repeatable manner. In the first stage, the algorithm calculates the AVR of a sample image taken at time t0 . For that, it locates the optic nerve and selects several radii of interest where the vessels are detected and measured. After that, the vessel segments are classified into veins and arteries and the AVR value is computed from a subset of the detected vessel segments [14]. The second stage allows to compute automatically the AVR in new sample images of the same patient’s eye acquired at times ti by registration. The image taken at time t0 is considered as the reference image and the AVR measurements in subsequent images are estimated in set of vessel segments selected for the reference image. The proposed methodology has been implemented as an AVR monitoring service on the SIRIUS system [15]. This system is a webbased application for retinal image analysis and patient’s checkup management. In this manner, we can assess the patient’s condition over time. The SIRIUS system also provides an automatic and a semiautomatic services to compute the AVR. At the semiautomatic mode, the system detects the vascular tree, measures the vessel widths and estimates the AVR. However, it still requires the user input to locate the optic nerve, select the vessels suitable for the calculus and classify them into arteries and veins. The paper is structured as follows. Section 2 explains the AVR monitoring method. In Section 3 the data set used is presented. Section 4 describes the performed experiments and shows the achieved results. Finally, in Section 5, we provide some concluding remarks. 2. Method The proposed method to estimate the AVR is based on a vessel registration approach which allow to measure the vessel widths at the same points from different sample images taken at different dates. This method consists of two main stages. In the first stage, a reference AVR is estimated using the application SIRIUS in a semior automatic way. This result is used to achieve subsequent AVR results by a vessel registration approach. The sample image taken at time t0 where the reference AVR is computed constitutes the reference image. In the second stage, we compute the AVR automatically in a new sample image from the same patient taken at time ti . First, this new image is registered to the reference image. Then, we locate in the new image the vessel segments in the reference image and we measure their widths. Finally, we apply the selected formula to get the AVR value. Fig. 1 shows a schema of the monitoring process. Fig. 2 includes two screenshots which show the integration of the AVR monitoring system in the SIRIUS web application. 2.1. Reference AVR measurement The reference AVR is measured by a module of the SIRIUS system [15]. The system has two modes of operation to estimate the
AVR. The first mode is semiautomatic and requires the user input to locate the optic disc, select the set of suitable vessels for the calculus and classify them into arteries and veins. The second mode is fully automatic, however, the user can supervise the result, discarding vessels, considering new ones or modifying the classification given by the system. In this section, we briefly describe the stages of the proposed automatic methodology. This comprises several stages. First, the optic nerve is located and then the vessels are segmented and measured in several circumferences concentric to the optic disc. After that, the vessel segments obtained are classified as arteries and veins. Finally, the system selects an appropriate set of vessel segments to estimate the AVR. 2.1.1. Optic disc location The optic disc location is based on the Hough transform [16]. First, we delimit the region of interest by a blob detection method, the Difference of Gaussians operator. For that, the original image is convolved with two Gaussian filters at different scales and the difference between these convolutions is obtained. The maxima of the DoG output corresponds to a rough location of the optic disc center. Then, we crop the image around this point using domain information and we apply the Hough transform. In order to restrict the set of points that can vote in the Hough transform, we detect the edges with the Sobel operator, whereas a crease detector algorithm is applied to obtain the vessel centerlines. If we consider the image as a topographic relief, where the pixel intensities represent the surface elevation, a crease is defined as the set of continuous image points that links the highest or lowest elevation levels. The crease detector used was the Multilocal Level Set Extrinsic Curvature enhanced by the Structure Tensor (MLSEC-ST) operator proposed in [17]. Thus, combining the Sobel and MLSEC-ST outputs, we obtain the optic disc perimeter points, these are, points that belong to edges but which are not near a vessel centerline. Once we have some optic disc perimeter points, we extract the optic nerve using the Hough transform for circle detection. 2.1.2. Measurement of vessel widths Once the optic disc is located, the vessels widths are measured in several circumferences centered on the optic disc. The number and radii of the circumferences can be changed by users in the semiautomatic mode. In the automatic mode, the circumferences are equally spaced a value of 0.5DD, being DD the optic disc diameter, and their radii range from 1DD to the border of the field-of-view (FOV). The proposed method to measure the vessel widths is based on snakes and it was originally proposed in [18]. A snake is a contour defined within the image which evolves guided by internal and external forces to fit some features. At first, the vessel centerlines are detected by the MLSEC-ST operator [17]. Thus, for a circumference of radius r, each snake consist of two chains of nodes initialized at the intersection points between each vessel centerline and the circumferences of radii r − n and r + n, with n = 10 pixels. Both chains of nodes advance along the normal to the vessel centerline in opposite directions toward the vessel walls guided by inner and outer energies. There are two types of nodes: common nodes which fit the edges and corner nodes which are adjusted to the corners of the vessel segment to be measured. In this model, the internal energy can be written as:
2 ∂s (s) ∂ s (s) Eint ((s)) = ˛(s) + ˇ(s) ∂s ∂ s2
(1)
where s (s) = [x(s), y(s)] defines the snake with node coordinates x(s), y(s), whereas the first and second term represent the first and second order derivatives, respectively. The parameters ˛(s) and ˇ(s)
S.G. Vázquez et al. / Computerized Medical Imaging and Graphics 37 (2013) 337–345
339
Fig. 1. Schema of the AVR monitoring process.
controls the snake shape. These values are constant for common nodes, in particular, ˛(s) = 0.25 and ˇ(s) = 0.01. The internal energy is null for corner nodes (˛(s) = ˇ(s) = 0) which implies the nodes are only attracted by the external forces. The external energy which attracts the snake to the target edge is defined by the next equation: Eext ((s)) = + ıEdist ((s)) + Egrad ((s)) + ωEstat ((s))
(2)
where is the dilation pressure, a vectorial magnitude which adjusts for each snake node its direction and sense of advance, Edist (i ) is the euclidean distance between a node and the vessel edges computed by the Canny operator from the input gray level image. These two terms are responsible for reaching the edges. On the contrary, the latter terms are stopping forces. Thus, the gradient energy, Egrad ((s)) brakes a node if its gradient is negative, that is, contrary to the advance direction. Finally, the stationary energy is based on the idea that probably a node should stop (keep on moving) if their adjacent nodes have stopped (kept on moving). This energy prevents movements of a node when its neighboring nodes have stopped or that a node stops in a local minimum when its neighboring nodes keep on moving. After the energy minimization, the snake contour should have a parallelogram shape and the nodes should be positioned in the vessel edges. However, image noise and discontinuities in the edge image may cause some nodes are placed outside or inside the edges. Hence an adjustment of snake contour is needed to correct the wrong nodes. Thus, a simple linear regression model is applied to the nodes at each snake side, achieving two regression lines. Then, the nodes are reorganized to fit to the corresponding regression line. If the correction is not possible, the snake is ruled out since it is inappropriate to be measured. In addition, if the final snake configuration has less than two nodes at each side corresponding to a vessel edge because many nodes evolved to the same position, the snake is also ruled out. Finally, the vessel width is estimated taken into account seven profiles in the snake parallelogram perpendicular to the centerline and equally spaced. At each profile, two distance measures are computed: the euclidean distance between the profile end points,
ω ˆ euc , and the distance obtained from a parabolic regression model, ω ˆ p . This last measure is used to correct the edge shift caused by a Gaussian smoothing applied to remove noise before applying the Canny edge detector. Thus, the profile gray levels are fitted to a parabola. To this end, only five points in the profile are considered: the end points, the point over the vessel centerline and the two points on either side of the latter with maximum gray level. Hence, the width of each profile is computed as the average of both measures (ω ˆ = (ω ˆ euc + ω ˆ p )/2). The final vessel width is estimated by the average of the three central values of the seven profile widths. In the automatic mode, before the vessel width measurement, the bifurcations and crossovers are detected [19] and they are ruled out. This step is necessary because both the measurement and classification of these points may introduce errors in the final AVR. 2.1.3. Vessel classification into artery and vein The classification of retinal vessels into arteries and veins is not straightforward. The biggest difference between arteries and veins is the color since veins are darker than arteries. However, the significance of this difference depends on the vessel size, being imperceptible for small vessels. In addition, the inter and intraimage lightness variability make the task more difficult. In [20], we have proposed an unsupervised method to classify the vessels which is robust to lightness and color variations. At each circumference, a coordinate axis centered on the optic disc divides the retina image into four quadrants. At each quadrant the vessel segments detected in the previous step are classified using the k-means clustering algorithm. Then, the coordinate axis is rotated through an angle of 20◦ and a new classification procedure is performed at each new quadrant. These two steps are repeated until it is achieved a 180◦ rotation. Then, the final class of a vessel segment is obtained, combining the results achieved in all the quadrants where the segment was classified. In order to avoid the influence of color variation, the input image is normalized using the multi-retinex technique [21]. Finally, we apply a procedure based on the minimal paths approach [22] to join vessel segments in consecutive circumferences which belong to the same vessel. Thus, the local results
340
S.G. Vázquez et al. / Computerized Medical Imaging and Graphics 37 (2013) 337–345
Fig. 2. Integration of the AVR monitoring system in the SIRIUS web application. (a) Interface for visualization and edition the AVR in the new sample image (first tab) and the reference image (second tab). (b) Interface which shows at right the list of the AVR measurements computed in new sample images using the reference AVR displayed at left side.
obtained in all the connected circumferences are combined to ensure the classification.
2.1.4. AVR computation Once the vessel segments have been classified and measured, we estimate the AVR as the ratio between the average artery and vein vessel widths [3,2]. Nevertheless, in the automatic mode, we make a vessel selection previously. In addition to the crossovers and bifurcations, we eliminate from the calculus the vessel segments whose width is over- or sub-estimated in relation with the rest of their connected vessel segments. According to experts, we do not take into account thin vessels, as well. Thus, we remove the arteries and veins whose widths are smaller than the 10-percentile [14].
Once the reference AVR has been estimated and saved for future computations, we also store some useful data for the vessel registration procedure, in particular, the optic disc center, the class of each vessel segment and the intersection points between each analysis circumference and the vessel centerline where the snakes have been initialized.
2.2. AVR measurement in a new sample image Given a new sample image taken at time ti , the method computes its AVR by registration using the stored reference AVR. For that, first, an image registration process between the two images is performed. Then the correspondence of the reference vessel
S.G. Vázquez et al. / Computerized Medical Imaging and Graphics 37 (2013) 337–345
341
Table 1 Image registration accuracy in a dataset of 20 retinal images.
Mean Standard deviation
a − r
MSE
RA(%)
0.033 0.050
0.097 0.070
99.100 1.550
and the number of seeds considered at each level were 6, 3 and 1, respectively. In the remaining levels, a Downhill Simplex algorithm is applied to optimize the correlation. At those levels, the search seeds are the optimized transformations of the previous level and the search process finishes when the difference between the maximum and the minimum values found in a neighborhood is lower than a preestablished threshold. Once the registration process has been finished and the two images are aligned, the registered crease images are used to obtain a similarity measure between them in order to determine the quality of the register process. This measure is the normalized cross-correlation coefficient, defined as:
Fig. 3. Schema of the multiple resolution registration process.
segments in the new sample image is obtained in order to measure their widths in the new image and compute the AVR in the same points. 2.2.1. Image registration The image registration algorithm used is a feature-based affine registration method whose landmark is the vessel tree [23]. Following the idea that the vessels can be thought as creases (ridges and valleys) when images are seen as landscapes, the registration process is based on the alignment of crease images. Thus, the first stage of the algorithm is to obtain the crease images of the reference and the new sample images using the MLSEC-ST operator [17]. After obtaining the crease images, an iterative optimization process to align both images is performed. That is, the reference image is fixed and the other one is transformed until a global maximum is accomplished. A suitable function to measure the alignment quality is the correlation function: Corr =
f (x)g((x))
(3)
x∈f
where f and g are the crease images and represent the transformation whose five parameters (x and y translation, tx , ty , rotation angle in clockwise direction and x and y scale, sx , sy ) we want to test. The search space is defined by the function Corr and the five parameters of the transformation. The function is non monotonic, that is, it has many local maxima and it is expensive to compute since it involves image transformations. Thus, to simplify the optimization process, the multiple resolution approach proposed by Elsen et al. [24] is carried out. Moreover, with the same aim in mind, the transformation is not applied to all pixels in the image, only the main creases with pixel values higher than a fixed threshold are transformed. The multiple resolution approach is handled by two pyramids (see Fig. 3) where the two crease images to register are at the bottom and each pyramid level is a half-resolution version of the image in the previous level, until images have a size of 64 pixels in each dimension. Then, the search starts at the top of the pyramid where an exhaustive search is accomplished and the correlation is computed in the Fourier domain. Thus, the search seeds are the image transformations which maximize Corr computed in the frequency domain. At this level, the transformations performed involve only rotations of an angle of 5◦ during a certain number of iterations. For a pyramid of three levels, the number of iterations was 10
=
x,y
x,y
[f (x, y) − f ] [g(x, y) − g ]
[f (x, y) − f ]
2
x,y
[g(x, y) − g ]
2
(4)
where g is the mean of the registered image and f is the mean of the reference image. The registration process is valid if is higher than a threshold, otherwise, the AVR cannot be computed using monitoring. In order to evaluate the registration method, we use a dataset of 20 retinal images taken with a Canon CR6-45NM non-mydriatic retinal camera with a 768 × 584 pixel resolution. For each image, 50 random transformations were applied with maximum translation and rotation values of ±100 pixels and ±5◦ , respectively, and without scaling. These transformations represent the maximum transformation values of the majority of the analyzed images. Table 1 shows the mean and standard deviation values of the absolute difference between the transformation applied and the transformation recovered by the method, the mean square error (MSE) of the recovered transformations and the registration accuracy (RA). The absolute difference was computed as a − r = |(tx − tx ) + (ty − ty ) + ( − ) + (sx − sx ) + (sy − sy )|, where a and r represent the applied and the recovered transformations being tx , ty , , sx , sy and tx , ty , , sx , sy their translation, rotation and scaling parameters, respectively. The registration accuracy represents the percentage of the transformations where the applied and the recovered transformations are equal with ±1 pixel of difference. To obtain this measure, four points are selected manually and their positions in the transformed and the dynamic images are compared. Thus, the average of the euclidean distances between the positions of the four points is calculated. If the absolute value of the average is less or equal than 1 pixel, the applied and the recovered transformation are aligned. The small MSE value and the high RA indicate the good performance of the image registration method. 2.2.2. Vessel registration Once a suitable transformation matrix between the reference image and new sample image has been obtained in the previous step and the transformation quality has been proved, we can measure the vessel widths in the new image at the same points considered in the reference image. First, the stored reference optic disc and intersection points are transformed into the new sample image (Fig. 4). For each transformed intersection point, we measure its distance, d, to the transformed optic disc center and we measure the width of the corresponding vessel segment by a snake. Once the snake evolution finishes, the topological adjustment is applied as well, and the vessel width is measured. The classification of the
342
S.G. Vázquez et al. / Computerized Medical Imaging and Graphics 37 (2013) 337–345
Fig. 4. Schema of the vessel registration procedure. In the new image (right), the vessel widths are measured at the intersection points registered from the reference image (left). Note that the same set of vessel segments is used to estimate the AVR.
vessels segments in the new image is obtained from the corresponding vessel segment in the reference image. Finally, the AVR is computed as the ratio between the average artery and vein vessel widths. Note that, we use the middle point of the snake initialization for the vessel registration and not the adjusted snake points. This is due to the fact that the final snake could be quite away from the initial position because of its deformation. If the snakes evolve from the same starting points in both images, there is a higher probability of similar deformations.
less than 2.2e−16 and a Pearson’s correlation coefficient of 0.61. Moreover, a p-value of 0.20 obtained in the second test reveals the AVR values are equal in average. However, the correlation coefficient and the p-value for the mean difference test increase to 0.92 and 0.89, respectively in a data set of 90 images, ruling out those images with the highest differences between the AVR values (see Table 2). The dispersion graphs of Fig. 5 show the correlation results. If we analyze the 25 images which present the AVR difference between expert and system highest (AVR6MMonitoring − AVR6MExpert ≥0.1) we derive the following assumptions:
3. Materials Eighty-nine hypertensive patients were studied from the Conxo Hospital in Santiago de Compostela (Spain). Two images, one per eye, have been acquired with a Cannon CR6-45NM non-mydriatic retinal camera for almost all patients. The images were stored with a resolution of 768 × 576 pixels in JPEG format. Each patient, after a basal assessment and an image acquisition, started an anti-hypertensive treatment. After 6 months, a new image was acquired. An ophthalmologist with many years of experience in the study of hypertensive retinopathy has analyzed the 158 pairs of images taken before and after treatment. The arteriovenous ratio in this set of images was computed by the expert using the SIRIUS system in a semiautomatic mode. In order to evaluate the monitoring system proposed in this work, each image taken at month 0 is considered the reference image and its AVR value graded by the expert represents the AVR reference. From the AVR reference, we compute the AVR in the image taken at month 6 and we compare the result obtained with the AVR graded by the expert in this image. 4. Results The proposed methodology has been tested in the data set composed of 158 pairs of images. We compute the AVR in the images taken at month 6 using the corresponding image captured at month 0 as reference image and its AVR value graded by the expert as AVR reference. A statistical analysis has been performed to test if the monitoring system behaves like an human expert, that is, whether despite using different vessels for the estimation, the AVR values obtained by the expert and the monitoring system are similar. We denote by AVR6MMonitoring the AVR obtained by the monitoring system and with AVR6MExpert the value measured by the medical expert in the images taken at month 6. A two-sided hypothesis test on Pearson’s correlation coefficient (H0 : AVR6MMonitoring ,AVR6MExpert = 0), as well as, a two-sided Welch’s t-test (H0 : AVR6MMonitoring − AVR6MExpert = 0) have been conducted. The first test provides evidence of correlation by means of a p-value
• The number of vessels taken into account to compute the AVR differs greatly from the image graded by the expert to the image measured by the monitoring system. This occurs mainly because the expert has selected a different number of vessels at months 0 and 6 as Fig. 6 shows or because the system does not detect some vessels at month 6 due to the image quality. As experts use a small number of vessels, these differences greatly affect the AVR value since it is computed as the ratio between the averages. • In other cases, it seems that the AVR value is different because the expert has selected in only one of the images a very characteristic vessel (wide or narrow) with a huge influence in the AVR (Fig. 6). • Moreover, an incorrect snake deformation is also possible, resulting in an under- or overestimation of the vessel width in one image whereas the width was properly measured in the other image. 4.1. Improvement of monitoring system First, we try to improve the system by preventing the elimination of vessels. The registration process works in the approved manner, so the vessel detection and measurement stages are the responsible of the elimination. Hence, we modify the measurement method, so each time a snake is invalidated due to its small number of nodes, it is thrown again, using as seed point a pixel located in an 8-neighborhood of the initial seed. The snake is definitely invalidated if any of the 8 possible seeds produces an enclosed snake. After introducing this change, the Pearson’s correlation coefficient in the 158 images increases slightly to 0.6329 while the p-value of the mean difference test remains the same, 0.209. Most of the missing vessels in the new sample image are lost because they are not detected. Their centerlines are not detected because the image quality. Thus, we modify the parameters of the centerline operator, which let to detect a greater number of vessels. These parameters are optimized to avoid false positives, but, since the vessel positions are known in the sample image, we can be less restrictive. Despite the change, the correlation value (0.6325) did not improve. This could happen because these parameters have
S.G. Vázquez et al. / Computerized Medical Imaging and Graphics 37 (2013) 337–345
343
Table 2 Test hypothesis results between the AVR values obtained by the human expert and the system in the images taken after 6 months of treatment. Pearson’s test p-value
Welch’s test p-value
0.616
