Optic nerve disc mathematical parameterization and the statistical comparison with measurements

The paper deals with the automatic measurement of the optic nerve disc diagnostic parameters from eye fundus digital images. The automation gives objective measurements, in contrast to measurements by physicians that have large variance. The disc shape has been detected by the use of mathematical morphological operations, of Canny’s edge detection, of geometric Hough’s transformation and of the approximation by the least square method. The statistical tests confirm that there is no significant difference between measurements performed by physicians and automatic measurements.


Introduction
Eye fundus examination is one of the most important diagnostic procedures in ophthalmology. A high quality colour photograph of an eye fundus helps in making diagnoses and follow-up of the development of a disease. Evaluation of eye fundus images is complicated and requires high-skilled experts for evaluation. The important parameter in diagnosis of glaucoma is an area of optic nerve disc (OND). However, as shown in this research, the measurements of OND, even performed by experts, have a rather large variance. Thus, objective parameterization of OND is crucial.
The parameterization was fulfilled by combination of the following image processing methods: mathematical morphological operations, Canny's edge detection, geometric Hough's transformation and the approximation using the least square method. The paper deals with some improvements to the former algorithm [9]. The statistical tests confirm that there is no significant difference between measurements performed by physicians and automatic measurements fulfilled by improved algorithm.

Localization of the optical nerve disc
There are three main structures in the eye fundus image: OND that appears in the eye fundus image as a yellowish disc, blood vessels and retina. The localization problem is extremely difficult because the blood vessels are located within the area of the OND and we do not even know a priori where the OND lies in the retinal image.
The retinal image is described as a set of three matrices I c = [I c xy ], c ∈ {r, g, b} for red, green, and blue components correspondingly, where I c xy is colour intensity at a pixel with coordinates (x, y). To achieve the homogenous structure (to remove the vessels), we have used the mathematical morphological closing operation that probe an image with a small shape (disc Z, described as a set of coordinates) known as a structuring element [5]. The diameter of Z (14 pixels) has been adjusted to the thickest vessel. The closing operation • is dilatation ⊕ followed by erosion Θ of brighter area: I c •Z = (I c ⊕Z)ΘZ, where I c ⊕Z = max (i,j)∈Z I c x+i,y+j , ∀x, y, I c ΘZ = min (i,j)∈Z I c x+i,y+j , ∀x, y. Successive image processing procedures are done individually with 2 matrices: 1) I p -the colour image recombined from closed one-colour images and transformed to grey scale [10]; 2) I g -the green colour image, as having little disruptive details. We will use the symbol I for both I p and I g , that are transformed individually parallel by the same methods. Also we will not change the symbol I after each stage of transformation. Further, the edges in the image were detected by the use of Canny's detector [1,2] consisting of these steps: • noise reduction by the use of the kernel smoothing where kernel is a twodimensional discrete Gaussian function; • determination of the edge points; • tracing edges to draw contours.
The points belonging to edges are determined analysing gradient norm ∇I xy = (D The angle of the gradient direction θ = arctan(D y I xy /D x I xy ) is rounded toθ nearest to nπ/4, n = 1, 2, 3, . . . . I is transformed by the use of the non-maximum suppression, id est, if a point is the local maximum of ∇I xy along gradient direction (comparing with two neighbours -alongθ and π +θ) it keeps its value I xy , otherwise I xy = 0.
With the tracing, the difficulty is that we cannot apply the fixed threshold level since there are no retinal images with identical properties. For automated threshold level calculation we use the Otsu's method [6,8]. The method calculates the optimum threshold separating two classes (bright and dark) so that their intra-class variances are minimal or, that is equivalent, the between-class variance is maximal. The calculated level is assumed as the upper level τ 1 and the lower level τ 2 = 0.1τ 1 . Two binary images T 1 and T 2 are produced applying calculated thresholds to the image after non-maximum suppression. We trace each segment in T 1 to its end and set them as contour points. At the segment end in the image T 1 we seek its neighbours in the image T 2 , since this image has much more details. If there are neighbouring pixels in the image T 2 , we denote them as contour points, too.
The next procedure is to locate OND among all contours. The circle shape with unknown centre (x c , y c ) and radius r is looked for by the use of Hough transformation [4]. The infinite number of circles may be drawn that cross one pixel; the centres of them constitute a circle around the pixel. If several such circles intersect at one point, then the point is the centre of the common circle looked for. Varying the r, the centre (x c , y c ) with the maximum number of approximate intersections and points  that constitute the contour of OND can be found. Finally, such points jointly from transformed gray-scaled I p and green band I g images are approximated by an ellipse according to the least square algorithm with constraints proposed by Fitzgibbon [2]. The OND ellipse and the original image are presented in Fig. 1.

Experimental tests
The optic nerve discs in s = 6 eye fundus digital images were measured by h = 15 experienced physicians in Kaunas University of Medicine (KUM) [7]. Images were presented on a computer monitor by an interactive programming tool and the physicians marked some points on the disc edge. The points were approximated by an ellipse to obtain parameters usual for physicians. Two physicians have not record their measurements of one of the images, so we have 88 measurements. The reorganisation and financial problems stopped experiments in KUM for a time therefore we may present here only initial results. The important parameter in diagnosis of glaucoma is an a -area of OND (area of the ellipse); more exactly, the decrement of the healthy OND area in time may indicate glaucoma. The measurements of OND area in pixels are presented in Fig. 2. It can be seen from the figure that for the same image the measurements by different physicians differ significantly, therefore, the objective measurements obtained using the algorithm is preferable. The measurements by the algorithm are presented in Fig. 2 and named as "Robot1". Pearson's correlation of measurements by a couple of physicians with other physicians' data are mostly less than 0.9, other physicians presented more correlated data, only one correlation coefficient of all matrix is less than 0.7. Robot1 data correlation coefficients with physicians' data are similar to the worst rows of physicians' data correlation matrix.
From the results presented in the Fig. 2 we see that measurements of OND area by Robot1 are shifted to large values in most cases. It is not a great disadvantage because the decrement of the area in the course of time is most important for diagnosis. However, for the first suspicion of the disease more exact measurements are important. The shift arises because physicians take into account not only yellowish disc but also bends of vessels inside and near the side of the disc. The edge detection is incapable to do that, considerably more sophisticated methods are necessary for that. We have proposed to decrease the shift by simple calibration of one parameter.
The relative errors of measurements were calculated by the formula r ij = (a ij − a j )/ā j , where i = 1, h is the index of the physician, i = 0 for Robot1; j = 1, s is the index of the image;ā j = h i=1 a ij /h is the mean of OND area measured by physicians -an "exact" area of jth OND. To decrease the shift of measurements, the axes of approximating ellipse, estimated by Robot1, are multiplied by the reduction coefficient k; so, the OND area is multiplied by k 2 . The sum of squared errors r 2 0j for j ∈ J ⊂ {1, . . . , s} is minimised seeking to adjust k: The solution of df (k)/dk = 0 gives us the minimiser k = a0j /āj a 2 0j /ā 2 j . The calibrated algorithm is named "Robot2". It is obvious that the multiplication of axes by k does not change the Robot2's data correlation with the physicians' data; however, it reduces the average of measurement errors. The averages of errors and corresponding values of k, adjusted for different subsets J, id est, adjusted to measurements by physicians for different subsets of images, are presented in Table 1. The k = 1 is for the Robot1. The question is can we distinguish the results of the algorithm from the measurements of physicians. Since we had only a few measurements, we have compared the distribution of all algorithm errors r 0j , j = 1, s with distribution of all physicians' errors r ij , i = 1, h, j = 1, s. The results of the significance tests are presented in the Table 1 were hypothesis H 0 represents some aspects of equivalence of the distributions.

Conclusions
The developed algorithm locates an optic nerve disc. The preliminary experiments indicate that the proposed calibration reduces the average of relative errors. The used statistical tests show that there is no significant difference between measurements obtained by physicians and by the calibrated algorithm; p values are greater than the level 0.05.