Quantification and characterisation of arteries in retinal images

Computer Methods and Programs in Biomedicine 63 (2000) 133 – 146 www.elsevier.com/locate/cmpb

Quantification and characterisation of arteries in retinal images Xiaohong W. Gao a,*, Anil Bharath b, Alice Stanton a, Alun Hughes a, Neil Chapman a, Simon Thom a a

Clinical Pharmacology, Imperial College School of Medicine at St. Mary’s, Paddington, London, W 2 1NY, UK b Centre for Biological & Medical Systems, Imperial College of Science, Technology & Medicine, London, UK Received 22 March 1999; received in revised form 9 March 2000; accepted 21 March 2000

Abstract A computerised system is presented for the automatic quantification of blood vessel topography in retinal images. This system utilises digital image processing techniques to provide more reliable and comprehensive information for the retinal vascular network. It applies strategies and algorithms for measuring vascular trees and includes methods for locating the centre of a bifurcation, detecting vessel branches, estimating vessel diameter, and calculating angular geometry at a bifurcation. The performance of the system is studied by comparison with manual measurements and by comparing measurements between red-free images and fluorescein images. In general an acceptable degree of accuracy and precision was seen for all measurements, although the system had difficulty dealing with very noisy images and small or especially tortuous blood vessels. © 2000 Elsevier Science Ireland Ltd. All rights reserved. Keywords: Retinal image; Vessel segmentation; Vessel geometry; Vessel diameter; Centre of the bifurcation; Image processing; Mathematical morphology; Gaussian function

1. Introduction Hypertension and associated cardiovascular risk factors cause structural alterations to blood vessels [1]. In general blood vessels are not visible in vivo, but the retina is an exception since, here, vessels can be seen through the pupil lying on the * Corresponding author. Present address: School of Computing Science, Middlesex University, Bounds Green Road, London, N11 2NQ, UK. Tel.: +44-181-3625000, ext. 2252; fax: +44-181-3626411. E-mail address: [email protected] (X.W. Gao).

surface of the retina, where they form a (essentially) two dimensional network. Viewing the retina is a common clinical procedure, which is useful in the diagnosis of hypertension, diabetes and several other diseases. Hypertension is associated with alterations in vessel diameters, branch angles and vessel length [2], and a number of studies have employed image analysis techniques to assess hypertensive retinal changes [2–4]. Some aspects of vessel topography can be described by the angle subtended by two daughter vessels at each bifurcation, or the junction exponent (x) which relates the diameters of vessels at a particu-

0169-2607/00/$ - see front matter © 2000 Elsevier Science Ireland Ltd. All rights reserved. PII: S0169-2607(00)00082-1

134

X.W. Gao et al. / Computer Methods and Programs in Biomedicine 63 (2000) 133–146

lar bifurcation using Murray’s law, namely d x1 + d x2 = d x0 (where d0, d1, and d2 are the diameters of parent and two daughter arterioles, respectively). These parameters may be useful descriptors of vascular changes as they influence functional properties such as power losses, shear stresses, and vessel density across a network. Accurate measurements are required to derive these parameters. Previously, measurements have been obtained manually by digitization of retinal images in combination with standard image analysis software JAVA (Jandel Scientific California). Such techniques are heavily user-dependent and time-consuming. Methods for more accurate measurement of vessel diameters have been described. Most involve improving the accuracy of edge location on an intensity distribution curve derived from a vessel cross-section [5]. Significant improvements were reported following the use of a single Gaussian function to model the vessel intensity profile [6]. However in high-resolution retinal images a single Gaussian function often is unsatisfactory to describe the intensity profile across a blood vessel, since in many cases a streak of light is clearly visible along the centre of the vessel. This is often termed the ‘central light reflex’ by clinicians and is believed to be due to reflection of incident light from the surface of the blood vessel wall and or the column of blood [7]. In this paper, we describe a computerised system for arterial quantification in retinal images. This includes vessel diameter estimation based on modelling of vessel intensity profile using twin Gaussian functions, vessel angle calculation and length measurements.

2. Materials For the purposes of the study, photographic negatives were obtained using a fundus camera (fx-50R, Kowa, Japan). Images were digitised from 35 mm negatives using a film scanner (LS1000, Nikon 35 mm film scanner) at a resolution of 2400¡ 2810 pixels. Two kinds of images were used: fluorescein and ‘red free’. Fluorescein angiography provides high quality images of the retinal vasculature. It yields high contrast between

the vessels and the background retinal layer, with sharp vessel edges allowing diameters to be determined readily and accurately. However, it is a relatively invasive technique involving intravenous injection of sodium fluorescein. ‘Red-free’ images do not involve the use of contrast agents, but a green filter is used to enhance the contrast between vessels and background. The image measurement system was mainly designed for use with ‘red-free’ images and was developed using C programming language with Xt/Motif library as a single process on a DEC Alpha workstation running UNIX operating system.

3. Methods and computerisation Fig. 1 shows the interface developed in the study and based on the MIDAS (Medical Information Display and Analysis System). It consists of two elements, one is a user-directed system and the other is the automatic measuring system. The user-directed system facilitates the manual measurement of angles, length and distance. For example, when a user clicks the button ‘Angle’ in Fig. 1, a message box pops up asking users to draw two angle lines using a mouse. Then the angle value in degrees will be displayed in the ‘Output’ window at the bottom-left graph of Fig. 1. Once the ‘MEASURE’ button is clicked, automatic measurement starts. Fig. 2 schematically illustrates the steps taken by the automatic system in measuring a vascular tree. It works on both red-free and fluorescein images with resolutions ranging from 1200¡ 1240 to 2100¡ 2400 pixels. The top-left quadrant of Fig. 2 is a red-free retinal image with a region of interest (ROI) selected by the user for vessel measurement. Users are then required to select two points using a mouse — the first close to the centre of the bifurcation of interest and the second close to the proximal bifurcation which gave rise to the parent vessel leading to the bifurcation of interest. The system then gives estimates of average diameters for each vessel (top-right quadrant), draws three lines along the centre of the vessels for calculation of the angular geometry at the bifurcation (bottom-left quadrant), and also

X.W. Gao et al. / Computer Methods and Programs in Biomedicine 63 (2000) 133–146

measures the length of the parent vessel from the proximal bifurcation to the bifurcation of interest (bottom-right quadrant). The strategy of this system involves refining the original estimate of the centre of the bifurcation provided by the user, detecting three vessel branches arising from the bifurcation, tracking each of these vessels to allow measurement of average diameters, calculation of angles between vessels and measuring the length of the parent vessel. The whole process takes about 1 min depending upon the length of the vessel to be measured. A flow chart illustrating

135

the strategy of the system is shown in Fig. 3 and is explained in detail below.

3.1. Refinement of bifurcation centre The precise location of the centre of the bifurcation is not only important for vessel tracking, but also for angle calculation: one pixel of difference in angle vertex position will result in two to three degrees of difference in calculated angle. The initial estimate provided by the user acts only a start point in the determination of the centre of

Fig. 1. The interface of MIDAS system.

136

X.W. Gao et al. / Computer Methods and Programs in Biomedicine 63 (2000) 133–146

Fig. 2. Image measuring system for quantification of retinal vascular network. The top-left quadrant shows a red-free image to be measured. The subsequent measurements for vessel diameter, vessel angle and vessel length are shown on the top-right, bottom-left, and bottom-right quadrants, respectively.

the bifurcation and small variation in this estimate is designed to have no effect on the ultimate location of the centre of the bifurcation. The following steps give the algorithm for refinement of the bifurcation centre. Centred at the given point, a region of interest (ROI) with size of 4d¡ 4d is created, where d is the initial estimation of the diameter for the parent vessel. It can be obtained either by modelling the intensity profile over the vessel cross section half way along the parent vessel as determined by the two points set by the user (the details of the modelling will be given later) or simply using a

fixed value, say 25 pixels, which represents the average diameter of the vessels in a group of similar images. The size of ROI is chosen to be large enough to enclose the whole bifurcation even if a users’ initial estimate is not close to the ‘true’ centre, but not too large to include other bifurcations. This region is then converted by optimal thresholding method [8] into a binary image using the following rules:

f(x, y)=

!

1 0

when f(x, y)] T otherwise

(1)

X.W. Gao et al. / Computer Methods and Programs in Biomedicine 63 (2000) 133–146

Where f(x, y) represents the ROI image and T the threshold. Fig. 4 demonstrates an example (left graph) of the binary image. Because of the central light reflex [7], there are some ‘holes’ with value of ‘0’ inside the vessel. Mathematical morphology of dilation (opening) followed by erosion transformations (closing) is therefore applied to fill those holes, i.e. for those pixels inside a vessel, replacing ‘0’ by ‘1’. The centre of the bifurcation is then found as follows: firstly, centred at each position of the binary image, a circle is drawn with all the pixels inside of the circle being ‘1’ and at least one pixel being ‘0’ outside of the circle. Secondly, the radius of this circle is calculated. The right graph of Fig. 4 is an example derived from the binary image of the left graph. The first estimate of the centre of the bifurcation is the location with the largest radius (r). If there are more than one such positions with the same value, the centre of those positions will be defined as the nominal centre of the bifurcation (BC). In the right graph of Fig. 4, the position with form of 19 is the centre of the bifurcation. In order to successfully define a unique centre for some noisy images, the above procedure for determination of BC is run twice using the first estimate of BC as the starting position for the second run. During testing in some very noisy images the algorithm visibly failed to correctly identify the location of BC. Consequently, after running the program a check was automatically performed by

137

calculating the distance between the user’s original estimate of BC and the new location. If the distance between these locations was greater than average vessel diameter (25 pixels in this group of images), then the user’s original estimate of BC was used and a default radius value of 15 was assigned for the subsequent processing stages.

3.2. Detection of 6essel branches The two positions initially chosen by the user only give information regarding the location of the parent vessel, the location of the two daughter branches also have to be determined. The strategy to achieve this is given below. Centred at BC (as determined above), a circle is drawn with a radius of 5r (r is the radius of the largest circle centred at BC, e.g. r is 9 (= 19/2) in the example of Fig. 4), which is illustrated in Fig. 5. When angle a is in a standard position and changes from 1 to 360° in the circle, the average intensity for all the pixels on each radial ray of the terminal side of the angle is calculated, together with the standard deviation (SD) of the intensities. For any radial ray lying within a vessel, the average intensity value on that line will be larger with a smaller SD for red-free images (with light vessel and dark background). For fluorescein images (dark vessels against light background), both average intensity and its corresponding SD will be smaller on the radius ray lying within a vessel area. Fig. 6 shows an example of average intensity values and their standard deviations (solid lines) plotted versus angle degrees for one complete cycle of 360°. This curve is then smoothed by a low-pass filter using a Sobel operator [8] along the x-axis using Eq. (2). fa=fa + 2+2fa+3+fa + 4−fa−4−2fa − 3−fa−1

(2)

The final curve (shown as dashed line in Fig. 6) is derived by further processing the above curve using the mean value of the SD and is given in Eq. (3). f(a) = Fig. 3. The flow-chart for the measuring system.

!

f(a) f(a−1)

0°5 aB 360°

if SD(a) B mean – SD otherwise (3)

138

X.W. Gao et al. / Computer Methods and Programs in Biomedicine 63 (2000) 133–146

Fig. 4. Binary images for a ROI transformed by an optimal threshold.

X.W. Gao et al. / Computer Methods and Programs in Biomedicine 63 (2000) 133–146

139

intensity calculated from the ray passing through such vessel is increased in addition to the average intensity. For some very small vessels with a3 B 5 pixels (a3 to be explained Section 3.3), it is not always possible to unequivocally define three peaks automatically. In these cases an observer’s intervention is required. When the system failed to locate three peaks, a dialogue appears asking the user to indicate two daughter branches using point and click with a mouse. The system then continues with the remaining tasks.

3.3. Vessel intensity profile modelling Fig. 5. The calculation of average intensity values on each angle radius line with radius of 5r.

where f(a) or fa is the average intensity curve processed using Eq. (2), mean – SD = SSD(a)/360, and f(0)=mean – SD. Three peaks can generally be identified which correspond to the three vessel branches. Occasionally there are adjacent vessels which pass through the ROI but do not arise from the bifurcation and so do not pass through BC. This approach successfully identifies these vessels as not belonging to the bifurcation since the SD of

Fig. 7 shows some examples of intensity distributions from different vessels of red-free images. In order to take account of the ‘central light reflex’, the intensity distribution over a cross-section of a blood vessel is modelled using twin Gaussian functions. One Gaussian function corresponds to the column of blood filling the vessel and the other corresponds to the central light reflex. The sum of the two functions therefore closely matches the vessel intensity profile. The intensity profile model is expressed in Eq. (4) I(x)=g1(x)− g2(x)

(4)

Fig. 6. Average intensity values and their standard deviations plotted with angles rotating from 0 to 360°. Curve with dashed line shows the final processed average intensity.

140

X.W. Gao et al. / Computer Methods and Programs in Biomedicine 63 (2000) 133–146

where x is the geometric distance and I(x) the intensity function. In practice, both g1(x) and g2(x) are assigned Gaussian function forms as given in Eqs. (5) and (6), respectively. 2

g1(x) =a1e − ((x − a2)/a3) +a4 g2(x) =a5e − ((x − a6)/a7)

2

(5) (6)

The parameter a1 of Eq. (5) gives the amplitude of the Gaussian function, a2 locates the peak position of the curve along x-axis, while a3 indicates the spread of the Gaussian curve and a4, the grey level of the background. The meaning of a5 – a7 is similar to the meaning of a1 – a3, respectively except for the central light reflex. The geometry of the function is representative of the geometry of vessels. For example, a2 represents the central position of the vessel and a3 is an estimator of the vessel diameter. For each vessel intensity curve derived from a normal cross section, parameters a1 – a7 in Eqs. (5) and (6) are calculated using a best-fit curve. Fitting is carried out by the application of non-linear Levenberg-Marquart method [9]. From a given set of proper initial values for the unknown parameters, this method works iteratively to minimise a x2 merit function and to determine the best-fit parameters. Fig. 8 gives two examples showing the performance of the vessel profile model. It shows that twin Gaussian functions provide more Fig. 8. Best-fit curves using twin Gaussian functions for vessel intensity distribution curves.

realistic estimates of the intensity profile across a vessel for this group of images.

3.4. Vessel tracking

Fig. 7. Intensity distributions over a cross section for different vessels.

Vessel diameters vary over different cross sections along a vessel. To obtain the average width, vessel tracking has to be performed for a given distance. Vessel tracking is also necessary to allow the vessel centre to be located at each cross-section to provide a centre line for angle calculation. In this study, vessel tracking consists of two steps. The first step is vessel segmentation and the second step is vessel measurement.

X.W. Gao et al. / Computer Methods and Programs in Biomedicine 63 (2000) 133–146

3.4.1. Vessel segmentation An arteriolar vessel is considered to be a linkage of many small vessel segments [10]. Each segment is identified by its centre location and its direction that is perpendicular to its cross sections. In other words, each small segment of a vessel is treated as a relatively straight cylinder. A vessel is then segmented according to the change in segment direction. In this study, each segment contained six successive searching points. Fig. 9 illustrates the process of vessel segmentation including allocation of segment centre (ci ) and determination of its direction (d( iri). The vessel segmentation process begins at BC (c0) with the vessel direction guiding point (cn ) obtained as described in Section 3.2 (or input by a user for the parent vessel). A normal line is drawn at a centre-search-point (CSP) ( in Fig. 9) with r distance away from BC (to avoid bifurcation area) and is perpendicular to the direction d( ir0 = (c0, cn ). The length of the normal line is fixed to 3r. Modelling of intensity distribution along the vector line using Eq. (4) then takes place. After modelling, CSP is re-positioned at the centre of Gaussian curve (a2 in Eq. (5) and shown as in Fig. 9) as one centre point of the vessel. The next CSP ( in Fig. 9) is then worked out based on current new centre point using Eq. (7): Distance (current-centre, next-CSP)= 2 pixels Direction (current-centre, next-CSP) = current vessel direction (7)

Fig. 9. Vessel segmentation process.

141

In Eq. (7), 2 pixels are the smallest integer sufficient to make the ‘next-CSP’ different from ‘current-centre’. The first part of vessel is considered to contain six successive centre points forming the first segment of the vessel. The centre of it (c1) is defined as the middle point of these six points. The d( ir1 = (c0, c1) is then taken as the direction of the first segment. Hence for segment i (i] 1), its direction is d( iri = (ci − 1, ci ), where ci is the centre of segments i. Therefore, the first segment has direction of (c0, c1) and the last segment has direction of (cn − 1, cn ). The direction of the previous segment is also the initial tracking direction for determination of the orientation of the next segment. A segment was defined as six successive centre points on the empirical basis that vessel direction does not change significantly within 12 pixels of distance (2 pixels by 6) and of that inaccuracies in any single centre point will not unduly affect the estimate of the segment centre (ci ). Vessel tracking is finished when it reaches the ending point cn or a specified distance of vessel has been tracked.

3.4.2. Vessel measurement After vessel segmentation, vessel measurement is performed by tracking every segment based on its new direction. As shown in Fig. 9, the initial direction (d( iri ) for the first segment represents an approximation based on the two ending positions of a vessel and does not represent true vessel direction. If the vessel direction is not correctly determined the normal line of the direction will not be perpendicular to the segment and the diameter estimate will be incorrect. A new direction is therefore given to the ith segment by d( iri (= (ci, ci + 1)) instead of (ci − 1, ci ) used in vessel segmentation. The last segment will be the (n− 1)th segment with the direction of d( irn − 1 = (cn − 1, cn ). Since cn is also a guiding position, tracking for vessel segmentation takes longer than vessel measurement. For example, during the measurement step, a vessel is tracked by a distance of 5r, while the distance for vessel segmentation takes about 7r.

142

X.W. Gao et al. / Computer Methods and Programs in Biomedicine 63 (2000) 133–146

Geometrically, this factor corresponds to the limits defining the area under the Gaussian curves as shown in Fig. 8 and can be used to locate vessel edges. In practise this value is fixed for a group of similar images. It can be estimated from comparisons with paired fluorescein images or manually determined diameters. Fluorescein images exhibit very sharp vessel edges, so that diameters obtained from this group of images probably provide a better estimate of vessel diameter values. Diameters can be measured effectively by applying a Sobel gradient to the vessel intensity profiles which exhibit the linkage of three straight lines instead of curved Gaussian shapes. The results of a comparison of paired red-free and fluorescein images are shown in Fig. 10. The slope of the linear regression (b) is the scaling factor and was calculated from the equation shown in Eq. (8). y=bx

(8)

From the data shown in Fig. 10, the scaling factor of b was calculated to be 2.47. This factor was used in subsequent determinations of diameters of red-free image as explained in Eq. (9). vessel-diameter=2.47 average (a3)

(9)

3.6. Angle calculation

Fig. 10. Comparison of measurements by the computerised system for both red-free and fluorescein images.

3.5. Vessel diameter determination A number of quantitative features of vessels including the location of vessel edges, diameter of a vessel, and vessel centre position are valuable for retinal vascular analysis [2,11]. Vessel diameter is estimated by multiplying a scaling factor to the mean vessel diameter estimators (a3 in Eq. (5)) obtained during vessel tracking for vessel measurement (second time of tracking).

At a bifurcation, the centre lines on two daughter branches and on the parent branch form an angular geometry [12]. As vessels may be tortuous (as discussed in Section 3.4), the decision where to draw a centred straight line for angle calculation is not straightforward. Experiments were therefore carried out to study how manual users measured angles. Most users tended to draw a vessel centre line along segments with similar directions. Therefore, in this system, a vessel centre line is formed by the best fit line derived using the first eight centre points within the first two segments plus BC. The best fit line is not constrained to pass BC, but 1/3 of the weighting is given to BC and 2/3 to the vessel points on the grounds that the crossing point of two vessel centre lines does not necessarily overlap BC, but is very close to it. Those centre points of a vessel are saved during

X.W. Gao et al. / Computer Methods and Programs in Biomedicine 63 (2000) 133–146

the procedure of vessel tracking for vessel measurement in Section 3.4. Angles are then calculated by applying the law of cosines of trigonometry shown in Eq. (10). c 2 = a 2 +b 2 − 2ab cos ÚC [ ÚC =cos − 1{(a 2 +b 2 −c 2)/2ab}

(10)

where a, b, and c are the length of each vessel centre line of a triangle and ÚC is the angle facing side c at vertex that is the cross point of two lines with length of a and b, respectively. Fig. 11 gives an example of an image illustrating the position of the three lines used for angle calculation.

3.7. Length calculation Vessel length is defined as a cumulative sum of the distance between two successive vessel segments as given in Eq. (11). Length =S distance(ci, ci + 1)

(11)

143

Table 1 Comparison of angle measurements (degrees) between the computerised system and manual methoda Bifurcation no.

Com

Man

Dif

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22

71.0 63.0 75.8 64.0 76.4 83.7 92.4 80.1 89.9 95.1 92.2 118.0 85.0 99.4 69.2 46.3 74.7 86.1 77.8 75.1 88.0 71.0

69.3 62.5 75.9 62.8 69.4 77.3 106.0 76.6 96.2 95.7 90.0 113.0 90.0 106.0 83.5 46.6 77.1 82.0 69.2 76.8 93.7 69.3

1.7 0.5 −0.1 1.2 7.0 6.4 −13.6 3.5 −6.3 −0.6 2.2 5.0 −5.0 −6.6 −14.3 −0.3 −2.4 4.1 8.6 −1.7 −5.7 1.7

80.6 14.5

81.3 15.9

−0.7 5.9

Mean SD a

In the table, ‘Com’ is the measurement by the computer system and ‘Man’ is manual measurement. ‘Dif’ is the difference between the two measurements.

where ci is the centre position of segment i, c0 is the BC and cn the ending point identified by the user. If the given ending point is at a bifurcation, it is refined to the centre of the bifurcation before the calculation as discussed in Section 3.3. This method not only takes account of the curvature of the vessel, but also avoids over-estimation of vessel length which would arise from calculating the distance between each individual centre point in each segment (there are six vessel centre points in each small vessel segment).

4. Results Fig. 11. An image showing the three vessel centre lines (black lines) used for angle calculation by the system.

Table 1 shows data comparing bifurcation angles measured manually or using the computerised

144

X.W. Gao et al. / Computer Methods and Programs in Biomedicine 63 (2000) 133–146

Table 2 Comparison of length measurement in pixels between the computerised system and manual methoda Bifurcation no.

L Com

L Man

L Dif

Cvt Com

Cvt Man

Cvt Dif

1 2 3 4 5 6 7 8 9 10

131.1 181.6 604.2 246.3 300.1 869.6 341.4 171.6 254.1 534.8

132.3 180.2 609 248.6 303.5 842.4 345.9 173.3 253.8 526.7

−1.2 1.4 −4.8 −2.3 −3.4 27.2 −4.5 −1.7 0.3 8.1

0.99 0.97 0.98 0.88 0.77 0.76 0.97 0.99 1 0.94

1 0.97 0.98 0.87 0.76 0.79 0.96 0.99 1 0.96

−0.01 0 0 0.01 0.01 −0.03 0.01 0 0 −0.02

Mean SD

363.48 235.1247

361.57 228.178

0.925 0.0864

0.928 0.084475

−0.003 0.012689

1.91 9.140509

a In the table, ‘Cvt’ represents curvature calculated by the length divided by the straight distance between two ending points, while ‘Com’, ‘Man’, ‘Dif’ have the same meaning as that in Table 1.

system. Manual data were calculated as the mean values of three separate measurements. The difference between the computerised measurements and manual was − 0.67° and within-subject standard deviation (sw) was 4.18° [13]. The comparison of length measurements obtained manually and by the computerised system is given in Table 2. The curvature of a vessel was defined as the length of a vessel divided by the straight distance between the two ending points of the vessel and is in the range of [0,1]. To illustrate the degree of curvature of a vessel, Fig. 12 shows a vessel tracked by the computerised system with curvature of 0.76. The black dots in the figure are the centre points of vessel segments and are used to calculate vessel length using Eq. (11). Manual measurement required observers to delineate the vessel by manually selecting points at intervals along the length of the vessel. The length was then calculated by adding the linear distance between each successive point. The positioning and interval between points was determined subjectively by the user. The difference between computer and manual determinations was 1.996.6 pixels (mean9sw) for length and −0.00390.009 (mean9sw) for curvature which indicates a good concordance between the measurements manually and those obtained using the computerised system.

Table 3 gives diameter measurements from a paired fluorescein and red-free image. The diameter of vessels in the fluorescein image was estimated using a Sobel operator as described above. The difference between fluorescein and red free measurements was 0.889 1.30 pixels (mean9 sw; n= 6). Unsurprisingly, larger deviations tended to occur for vessels with diameters B 10 pixels (a3 B 4).

5. Discussion Hypertension and diabetes are major causes of global morbidity and mortality the effects of which are manifest in the retinal circulation. Qualitative assessment of the retinal vasculature provides useful clinical information about disease progress and severity. However the information obtained from retinal examination is likely to be further improved by applying a more quantitative approach to the assessment of retinal vascular changes. One of the difficulties with such an approach is the subjective and time-consuming nature of many such measurements. This paper describes a computer-based semi-automatic system for the quantification of arterial branch parameters that are affected by hypertension [2]. In general an acceptable degree of accuracy and

X.W. Gao et al. / Computer Methods and Programs in Biomedicine 63 (2000) 133–146

145

precision was seen for all measurements, although the system has occasional difficulties dealing with noisy images and small blood vessels. Although some improvements can be anticipated from improved image acquisition either using improved film or direct digital image capture, determination of bifurcation angles remains difficult in some instances. This is particularly the case when vessels are very tortuous as is illustrated by Fig. 13 where the system fails to determine vessel centre lines as a result of the abrupt change in vessel direction and the validity of the angle measurement is dubious. Further work is needed to satisfactorily solve the problems posed by such vessels.

Fig. 12. An image showing the points used for vessel length calculation on a blood vessel with curvature of 0.76.

Table 3 Comparison of diameter measurements in pixels for arterioles from a group of fluorescein and red-free image pairs Bifurcation no.

Com

Man

Dif

1 2 3 4 5 6

39.5 34 18.5 23.3 22.7 9.1

39.1 34.5 18.5 23.9 22.9 13.5

0.4 −0.5 0 −0.6 −0.2 −4.4

Mean SD

24.5 9.94

25.4 8.8

−0.9 1.6

Fig. 13. System fails to draw vessel centre lines (right hand side) for vessels with direction changed abruptly.

146

X.W. Gao et al. / Computer Methods and Programs in Biomedicine 63 (2000) 133–146

6. Summary A computerised system has been developed to allow semi-automatic quantification of vessel geometry in retinal images. This system utilises digital image processing techniques to rapidly provide reliable and more comprehensive information for the retinal vascular network than can easily be obtained by manual measurements. Several strategies and algorithms have been developed including those for determination of the location of the centre of a bifurcation, detection of vessel branches, vessel tracking, estimation of vessel diameter and length and calculation of angular geometry at a bifurcation. Comparison of computer determined findings from ‘red-free’ images with manual measurements or measurement in fluorescein images indicates that, in general, the system gives satisfactory results in terms of accuracy and precision.

[3]

[4]

[5]

[6]

[7]

.

[8] [9]

Acknowledgements This project was funded by Medical Research Council (MRC), UK. Their support is gratefully acknowledged.

[10]

[11]

References [12] [1] A.D. Hughes, M. Schachter, Hypertension and blood vessels, Br. Med. Bull. 50 (1994) 356–370. [2] A.V. Stanton, A.B. Wasan, A. Cerutti, S. Ford, R. Marsh, P.P. Sever, S.A. Thom, A.D. Hughes, Vascular network

[13]

changes in the retinal with age and hypertension, J. Hypertens 13 (1995) 17 – 24. B. Khoobehi, G.A. Peyman, K.D. Vo, Relationship between blood velocity and retinal vessel diameter, ARVO Abstract, Invest. Opthalmol. Vis. Sci. 33 (4) (1992) 804. H.C. Chen, J. Wiek, S.M.B. Rassam, V. Patel, E.M. Kohner, Retinal vessel diameter changes in the cardiac cycle, ARVO Abstract, Invest. Opthalmol. Vis. Sci. 33 (4) (1992) 811. V. Patel, S.M.B. Rassam, O. Brinchmann-Hansen, O. Engvold, E.M. Kohner, A new concept in accurately measuring retinal vessel diameter from transmittance and densitometry profile of fundus photographs, ARVO Abstract, Invest. Opthalmol. Vis. Sci. 33 (4) (1992) 804. L. Zhuo, M.S. Rzeszotarsk, L.J. Singerman, J.M. Chokreff, The detection and quantification of retinopathy using digital angiograms, IEEE Trans. Med. Imaging 13 (4) (1994) 619 – 626. O. Brinchmann-Hansen, O. Engvold, Microphotometry of the blood column and the light streak on retinal vessels in fundus photographs, Acta Ophthalmol. 64 (supplement 179) (1986) 9 – 19. R.C. Gonzalez, P. Wintz, Digital Image Processing, Addison-Wesley, Reading, UK, 1987. W.H. Press, B.P. Flanney, S.A. Tenkolsky, W.T. Vetterling, Numerical Recipes in C — The Art of Scientific Computing, Cambridge University Press, Cambridge, UK, 1990. Y. Sun, Knowledge-based segmentation and correspondence of vascular structures from biplane angiograms, medical imaging IV: image processing, SPIE 1233 (1990) 257 – 265. American Academy of Ophthalmology, Ophthalmic Pathology, in: Basic and Clinical Science Courses, section 11, 1991, 179. A.G. Roy, M.J. Woldenbery, A generalisation of the optimal models of arterial branching, Bull. Math. Biol. 44 (1982) 349 – 360. D.G. Altman, Practical Statistics for Medical Research, Chapman & Hall, London, UK, 1991.