- Email: [email protected]
Proposal, verification and comparison of three computer image analysis methods for detection and evaluation of colour glaucomatous changes within the optic disc of a human eye retina
ARTICLE IN PRESS
Optik
Optics
Optik 116 (2005) 133–143 www.elsevier.de/ijleo
Proposal, verification and comparison of three computer image analysis methods for detection and evaluation of colour glaucomatous changes within the optic disc of a human eye retina Frantisˇ ek Pluha´cˇeka,, Jaroslav Pospı´ sˇ ilb a
Department of Optics of Palacky` University, 17. Listopadu 50, 77200 Olomouc, Czech Republic Department of Experimental Physics and Joint Laboratory of Optics of Palacky` University and Institute of Physics of Academy of Sciences, 17. Listopadu 50a, 77200 Olomouc, Czech Republic
b
Received 20 August 2004; accepted 28 November 2004
Abstract The typical symptom of the human eye glaucoma is a rise and a progression of the bright area (named pallor area) within the retina blind spot. The image analysis manner proposed by the authors detects and suitably numerically describes the relative size of the representative pallor area in the colour digital image of the retina obtained by a suitable fundus camera connected with the computer. Three new different computer image analysis statistical methods for experimental diagnostic evaluation of the obtained characteristic data are proposed in this article: the quantile curves method, the neural net method and the probability density curves method. The quantile curves method is based on the graphical comparison of a relative representative pallor area size with its determined normal value. The neural net and probability density curves methods can automatically and objectively classify the investigated eyes in exactly defined glaucoma risk classes and diagnosed glaucoma with the rated probabilities of incorrect diagnosis determination. All mentioned methods are verified and mutually compared by their application to the large statistical sets of human retina images of various healthy and glaucomatous subjects. r 2005 Elsevier GmbH. All rights reserved. Keywords: Glaucoma; Optic disc; Image analysis; Statistical evaluation; Neural net
1. Introduction Glaucoma is the common name of a group of single human eye diseases characterized by an increased intraocular pressure (the eyewater pressure in the eyeball), typical damage to the photosensitive cell nerve fibres of the retina, colour and shape changes of the blind spot of the retina (named optic disc, OD) and in Corresponding author. Tel.: +420 585634310.
E-mail addresses: [email protected] (F. Pluha´cˇek), [email protected] (J. Pospı´ sˇ il). 0030-4026/$ - see front matter r 2005 Elsevier GmbH. All rights reserved. doi:10.1016/j.ijleo.2004.11.009
the eye visual field [1–6]. The prevalence of glaucoma is approximately in 1–2% of adults [1]. The patient with this disease feels no subjective problems for a long time, but there is a slow decrease of the visual functions. The final stage of non-medicate glaucoma is practical blindness. Thus, the early diagnostics of glaucoma is very important. For the progression of this disease, the size of an iridocorneal angle can be representative, see publications [1–4]. The considered iridocorneal angle is the anterior chamber space between the iris base and the inner cornea margin. Its size influences the intraocular pressure and the glaucoma risk. In view of the discussed
ARTICLE IN PRESS 134
F. Pluha´cˇek, J. Pospı´ sˇ il / Optik 116 (2005) 133–143
angle size, two glaucoma types are commonly distinguished: the open-angle glaucoma (if the angle is of normal size) and the angle-closure glaucoma (if the angle is close). Only the routine open-angle glaucoma is considered in this article. One of the first symptoms of glaucoma is the decay of the retina photosensitive cell nerve fibres and a subsequent progress of the depression or cup within the OD. These glaucomatous processes cause a rise and enlargement of the bright area (named pallor area) within the OD cup bottom, which represents the OD area with the absence of the nerve fibres [1]. In the case of a healthy eye, the pallor area is missing or is only small. Because the subjective determination and evaluation of the pallor area or other symptoms have small reproducibility, the objective diagnostic methods are required. The size of the cup can be objectively measured and evaluated at present by the method of the computer evaluation of stereoimages [7–11] and by the laser scanning tomography method [12–15]. It is also possible to scan the human eye retina using a suitable colour fundus camera connected to the computer by a convenient graphic card, as was carried out by the authors of this article. A colour digital image obtained in this way can be advantageously archived in an appropriate data medium and can be used for detection of the colour glaucomatous changes using a suitable computer (digital) image analysis. For diagnostic purposes, the relative size of the so-called representative pallor area in the image must be considered because of its dependence on the considered so-called representative OD size. It is the ratio of the representative pallor area size and the representative OD area size in the image (denoted as P/D area ratio or only P/D) which is chosen as the relative pallor area size measure in the case of the utilized image analysis manner presented in this article. The representative OD area is conventionally divided into several suitable parametric sectors and the P/D area ratio is computed for each sector. For diagnostic evaluation of the P/D area ratios, three new statistical methods of computer image analysis are proposed, verified and compared in this paper. In accordance with the utilized principle and form of the obtained main results, the following names for these methods were accepted: the quantile curves method, the neural net method, and the probability density curves method. All computer calculations were realized by the mathematical programs MATLAB and MATHCAD.
2. Used equipment and the measuring conditions The colour computer digital OD images of the size 576 768 pixels, obtained by the computer-connected
digital fundus camera DFK 98 [16], were used. The mentioned camera contains the colour CCD videocamera connected with the adapted photographic system Retinofot Carl Zeiss Jena 201, see [16]. During the photographic scanning of the OD, the described equipment was located so that the rows of the acquired image matrix correspond to the horizontal direction. The actual size of every OD area was approximately 5.8 mm 4.4 mm2. The brightness values of each of the corresponding three image colour channels (i.e. red, green and blue channel) were quantized to 256 discrete levels and represented by integers from 0 to 255. Value 0 represents the minimum brightness and value 255 the maximum brightness of each colour digital image. All the following considerations and results relate to these measuring conditions. The maximum number of 240 OD images of 120 normal healthy subjects and of 240 OD images of 120 subjects with the open-angle glaucoma were utilized in our research. According to an eventual medical requirement, additional medical diagnoses were proposed by the cooperate glaucomaspecialized ophthalmologist on the basis of independent measurement of the intraocular tension, visual field, enlargement of the cup using the laser scanning tomograph and of a subjective observation of the thickness of the retina nerve fibre layer. The examined normal healthy subjects had the age range from 14 to 76 years with the average age of 46 years and with the standard deviation of 13 years. The glaucomatous subjects under examination had the age range of 10–82 years with the average of 51 years and the standard deviation of 16 years.
3. Performed preliminary determination of the representative OD area and the relative pallor area size by the utilized computer image analysis manner 3.1. Brightness of the representative OD area and definition of the representative pallor area The representative pallor area is situated within the representative OD area in the retina image obtained by the above-mentioned fundus camera. The representative OD area was manually demarcated by a convenient little excentric representative boundary ellipse, which connects the five auxiliary points selected by the operator’s experience in the used colour digitized image of the OD. In order to meet the requirement of objective orientation of the image, a polar coordinate system was chosen with its pole 0 in the centre of the representative boundary ellipse and with the polar line oriented horizontally according to the reverse direction of the human nose, see
ARTICLE IN PRESS F. Pluha´cˇek, J. Pospı´ sˇ il / Optik 116 (2005) 133–143
Fig. 1. The angles were considered clockwise for the right eye and anticlockwise for the left eye. It was ascertained experimentally that the OD image brightness G is most expressive in the green colour channel and is displayed with sufficient high optical contrast. Hence, only this channel was used for the needed normalized brightness determination. For such normalization the two suitable reference pixel groups A and B were considered in accordance with Fig. 1. It is evident that the image point (u,v) located in the intersection of the uth row and the vth column of the pixel array lies within these pixel groups of empirically chosen width a/5 when hold the inequalities qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a rðfÞ p ðu u0 Þ2 þ ðv v0 Þ2 prðfÞ; (1) 5 where 15 pfp 15 in the case of the group A or 165 pfp195 in the case of the group B. The point ðu0 ; v0 Þ represents the position of the pole 0 of the abovedefined polar coordinate system, a is the reference major semi-axis size of the representative boundary ellipse, f represents the oriented angle between the position vector of a pixel (u,v) and the polar line, rðfÞ is the mutual distance between the pole 0 and the boundary ellipse point of the angular coordinate f: Because the brightness values of the OD image areas corresponding to the pixel groups A, B are not practically influenced by the glaucoma progression, they can be used to determine the reference brightness values GA and GB for the already mentioned normalization of G. These values are defined as the green image channel brightness values of the largest occurrence in pixels of the groups A and B. Moreover, the location of each pixel group can be represented by the average values u¯ A ; v¯ A and u¯ B ; v¯ B : The also required normalized brightness Gn(u,v) in each point (u,v) of the representative OD image is proposed by the useful relation G n ðu; vÞ ¼
2Gðu; vÞ : GA þ GB
(2)
135
Fig. 2. Example of the healthy (left) and glaucomatous (right) OD image with the corresponding negligible small and expressive large bright pallor area. The boundary ellipses define the representative OD areas.
The above-mentioned values GA and GB differ in general, because the brightness of the OD periphery often inconveniently changes in the polar line direction. For the effective compensation of this effect, we designed and experimentally verified the acceptable brightness transformation G 0n ðu; vÞ ¼ G n ðu; vÞ G An
G Bn G A n ðv v¯ A Þ v¯ B v¯ A
(3)
of the normalized brightness Gn(u,v) in the pixel (u,v) to the new value G0n ðu; vÞ: The introduced values G An and G Bn are the values GA and GB, which are normalized by value (2). The corresponding representative pallor area of the OD image is then defined as the image area within the representative boundary ellipse, where the normalized transformed value G 0n is greater than or equal to the selected suitable threshold value p. The empirically established equality p ¼ 0:2 is accepted in the following text. Fig. 2 shows an example of healthy (left) and glaucomatous (right) OD image with the corresponding negligible small and expressive large bright pallor area. The boundary ellipses define the representative OD areas.
3.2. Relative pallor area size characteristic vectors and the standard deviations of their measurement
Fig. 1. Example of an optic disc representative boundary ellipse, the considered polar coordinate system and the chosen pixel groups A and B.
Let us assume the knowledge of the representative OD and pallor areas which were determined in agreement with the previous text. The relevant relative pallor area size P/D, mentioned in Chapter 1 of this article and related to the total representative OD area, can be conventionally numerically represented by the 36-dimensional characteristic vector (P/D1, P/D2,y, P/D36) or by the five-dimensional characteristic vector (P/DT, P/DI, P/DII, P/DIII, P/DIV) in accordance with [14,15]. The introduced single components P/D1, P/D2,y,P/D36 represent the elementary P/D area ratios for 36 equally spaced OD sectors with the common vertex in the pole 0 of the used polar coordinate system
ARTICLE IN PRESS 136
F. Pluha´cˇek, J. Pospı´ sˇ il / Optik 116 (2005) 133–143
and with the vertex angle ranges h0 ; 10 Þ; h10 ; 20 Þ; . . . ; h350 ; 360 Þ: On the other hand, the single components of the five-dimensional characteristic vector (P/DT, P/DI, P/DII, P/DIII, P/DIV) denote successively the P/D area ratios of the total representative OD area and of its four quadrants numbered I, II, III, IV and contained the vertex angle ranges h45 ; 45 Þ; h45 ; 135 Þ; h135 ; 225 Þ; and h225 ; 315 Þ: The performed experiments showed that the used subjective manual demarcation of the representative OD elliptical border for a given OD image is markedly influenced by a significant variability, which can be represented by the single standard deviations of the above defined and measured characteristic components. For the performed estimation of such measurement single deviations, the demarcation of the every representative OD elliptical border and also the subsequent computation of all the corresponding partial P/D area ratios were always repeated ten times for randomly selected 40 healthy and 40 glaucomatous OD images (i.e. 800 demarcations were performed). Then for each selected image, the obtained partial P/D area ratios were statistically processed in the way that the associated single standard deviation ss ; s 2 fT; I; II; III; IV; 1; . . . ; 36g; of the repeated measurement of a partial P/Ds area ratio is represented by the mean value of all corresponding partial P/Ds standard deviations. Furthermore, for practical usage, the single values s1 ; s2 ; . . . ; s36 can be often suitably replaced by their average of the possible expression s1;2;...;36 : Our computations yield the values sT ¼ 0:008; sI ¼ 0:02; sII ¼ 0:01; sIII ¼ 0:007; sIV ¼ 0:01; s1;2;...;36 ¼ 0:02: The results of investigated dependences of the P/DT discrete values on the age of the examined 240 healthy and 240 glaucomatous subjects were suitably approximated by the corresponding regression lines using the mathematical theory of regression analysis [17,18]. The appropriate regression lines of equations P=DT ¼ 5:6 105 age þ 0:1; P=DT ¼ 1:7 103 age þ 0:3 with the correlation coefficients 0.01, 0.3 and with the standard deviation 0.3 for both regression line slopes were obtained by computation. The used variable quantity under indication age means the subject age in years. The obtained small values of the regression line slopes and the relative high values of their standard deviations indicate that the age dependences are insignificant for the considered subject sets.
4. Quantile curves method for diagnostic evaluation of the 36-dimensional P/D area ratios The above-defined 36-dimensional characteristic vector (P/D1, P/D2,y,P/D36) of the P/D area ratio, which
is fundamental for this method, was computed for each of the 240 healthy and 240 glaucomatous OD images mentioned in Chapter 2. The obtained statistical probability distributions of the values of each single component P/Dn, n ¼ 1; 2; . . . ; 36; can be advantageously mathematically compared by the associated experimental distribution functions Dhn(P/Dn) in the case of healthy OD images and Dgn(P/Dn) in the case of glaucomatous OD images. According to [14], the special ð90Þ ð50Þ ð90Þ parametric values qð50Þ of P/Dn, hn ; qhn and qgn ; qgn fulfilling the general relations Dhn ðqðxÞ hn Þ ¼ x=100;
x 2 f50; 90g;
(4)
Dgn ðqðxÞ gn Þ ¼ x=100;
x 2 f50; 90g
(5)
and named as x% normal quantiles and x% glaucomatous quantiles, were utilized. Thus, the probability of occurrence of the P/Dn values in normal healthy eyes, or glaucomatous eyes, up to the extreme value ðxÞ qðxÞ hn ; respectively, qgn is x%: For our comparison, the ðxÞ ðxÞ ðxÞ ðxÞ ðxÞ considered values qðxÞ h1 ; qh2 ; . . . ; qh36 and qg1 ; qg2 ; . . . ; qg36 of P/Dn were at first well-arranged in the descending order and then plotted in dependence on their sequence in this ordering (see Fig. 3). The obtained graphs represent the x% normal quantile curves (solid lines) and the glaucomatous quantile curves (dashed lines) for the considered sets of eyes. Moreover, the 36-dimensional characteristic vector can be also applied only to one arbitrary eye. Such eye is graphically described by a one-eye curve which represents the dependence of its single P/Dn values sorted in the descending order on their sequence in this ordering. Examples of the
Fig. 3. The obtained 50%- and 90%- normal quantile curves (solid lines) in comparison with the 50%- and 90%glaucomatous quantile curves (dashed lines) and with oneeye curves of one randomly selected healthy (normal) eye (dotted line) and glaucomatous eye (dot–dashed line). These curves relate to the described quantile curves method.
ARTICLE IN PRESS F. Pluha´cˇek, J. Pospı´ sˇ il / Optik 116 (2005) 133–143
137
obtained one-eye curves for one randomly selected healthy (normal) eye (dotted line) and glaucomatous eye (dot–dashed line) are also shown in Fig. 3. It is evident that the values of P/D area ratios, presented in the graphical form for the glaucomatous eye, are markedly influenced by the glaucoma. Generally, if the one-eye curve of an investigated eye is above the 90% normal quantile curve, the high risk of glaucoma occurrence in this eye can be assumed. Such a graphical comparison can be a useful glaucoma diagnostic tool.
5. Neural net method for diagnostic evaluation of the five-dimensional P/D area ratios 5.1. Description of the structure of the used artificial computer-simulated neural net Generally, the artificial neural net represents the mathematical transformation system with a N-dimensional vector x ¼ ðx1 ; . . . ; xN Þ of input components and a M-dimensional vector y ¼ ðy1 ; . . . ; yM Þ of output components, see publications [19,20]. The basic element of every neural net is a formal neuron, whose inputs contain the information from other neurons or from the outer environment. Each neuron has only one output, transmitting the resulting information from this neuron to further elements (neurons) or to the output of the net. All the elements are arranged in L layers and each neuron in the lower layer is connected with each neuron in the layer above it. The neurons in the same layer are not mutually connected. Inputs of the neurons of the lowest layer (denoted as the first or input layer) are real inputs of the net and outputs of the neurons of the highest layer (denoted as the Lth or output layer) represent the real outputs of the net. According to Refs. [19,20], the mathematical formal model of the kth neuron in the ith layer can be represented by the equation " !#1 N i1 X xk;i ¼ 1 þ exp wi1;j;k xj;i1 yk;i ; (6) j¼1
where xk,i is the output of the kth neuron in the ith layer, xk;0 ¼ xk is the kth input of the whole neural net, xk;L ¼ yk is the kth output of the whole neural net, wi,j,k represents the weight of connection of the jth neuron in the ith layer with the kth neuron in the (i+1)th layer and w0,j,k is the weight of connection of the jth neural net input with the kth neuron in the first layer ði ¼ 1Þ: Furthermore, Ni denotes the number of neurons in the ith layer, N 0 ¼ N is the number of the neural net inputs (called the dimension of the input), N L ¼ M is the number of the neurons in the output layer (called the dimension of the output) and yk;i represents the
Fig. 4. Scheme of the considered neural net structure for diagnostic evaluation of the five-dimensional P/D area ratios.
threshold of the kth neuron in the ith layer. It comes true that the output signal of each neuron in the ith layer represents the input of all neurons of the (i+1)th layer at the same time. The neural net with the parameters L ¼ 2; N ¼ N 0 ¼ 5; N 1 ¼ 3; M ¼ N 2 ¼ 1 and with the input vector x ¼ ðP=DT ; P=DI ; P=DII ; P=DIII; P=DIV Þ is considered in the following text. Its structural scheme is shown in the Fig. 4.
5.2. Back-propagation learning algorithm The values of the already mentioned thresholds and weights of the considered neurons in the neural net were determined by the so-called back-propagation learning algorithm, see [19,20]. In accordance with this algorithm, the thresholds and weights are iteratively adapted to minimize the error function E xð1Þ ;xð2Þ ;...;xðCÞ ;d ð1Þ ;d ð2Þ ;...;d ðCÞ ðw; hÞ ¼
C X M 1 X 2 ðyðiÞ d ðiÞ j Þ ; 2 i¼1 j¼1 j
(7) where w is the vector containing all weights in the neural net and h is the vector containing all thresholds of all ðiÞ neurons. The quantities yðiÞ ¼ ðyðiÞ 1 ; . . . ; yM Þ; i ¼ 1; . . . ; C; are the real output vectors of the neural net responding ðiÞ to the given set of the input vectors xðiÞ ¼ ðxðiÞ 1 ; . . . ; xN Þ ðiÞ (denoted as a learning set), d ðiÞ ¼ ðd 1ðiÞ ; . . . ; d M Þ are the required (ideal) output vectors corresponding to the input vectors xðiÞ and C represents the range (number) of the learning set. The chosen learning set contains 240 input vectors constructed using the P/D area ratios of 120 selected healthy OD images and of 120 selected glaucomatous OD images with the required neural net output binary values equal to the minimal, respectively,
ARTICLE IN PRESS 138
F. Pluha´cˇek, J. Pospı´ sˇ il / Optik 116 (2005) 133–143
maximal possible value 0, respectively, 1 of the neural net output. The considered iterative learning algorithm starts with the initial values w(0), hð0Þ: Simultaneously, the single components of the initial vectors w and h should be low [19,20], e.g. chosen at random from the interval (0.3, +0.3). In the lth step of the considered learning process, the values of weights and thresholds change from the actual values w(l) and hðlÞ to the new values w(l+1) and hðl þ 1Þ according to the relations wðl þ 1Þ ¼ wðlÞ Z1 gw þ a1 ½wðlÞ wðl 1Þ;
(8)
hðl þ 1Þ ¼ hðlÞ Z2 gy þ a2 ½hðlÞ hðl 1Þ;
(9)
where the symbols gw, respectively, gy represent the vectors of partial differentiations of the error function (7) with respect to the weights, respectively, thresholds, of the neurons in all layers of the neural net and Z1 ; Z2 ; a1 ; a2 are suitable learning parameters. We realized 50 learning iterations with empirically determined values Z1 ¼ Z2 ¼ 0:1; a1 ¼ a2 ¼ 0:005: The dependence of the error function (7) on the number of iterations of the utilized learning process is shown in Fig. 5. Simultaneously, the used values of w and h; also named the learned values, relate to the considered last iterative step. After the above-mentioned learning, the considered neural net was tested [19,20] using the test set of 240 input vectors constructed from the P/D area ratios of 120 selected healthy eye images, 120 selected glaucomatous eye images and the corresponding required output values that were not used in the learning set. According to the considered learning set mentioned above, the required neural net output was chosen equal to 0 for healthy eye images and 1 for glaucomatous eye images. During the test, the value of the error function (7) was computed for both the introduced test set and the used learning set. The obtained difference DE between the
Fig. 5. Dependence of the error function (7) on the number of iterations of the utilized learning process.
error function value for the test set and the learning set was acceptably small, DE 0:18: In view of the nonzero single standard deviations sT ; sI ; sII ; sIII and sIV ; the corresponding mean standard deviation sn 0:02 of the neural net outputs was numerically estimated on the basis of the common error-propagation rule [17].
5.3. Proposal of the neural net eye diagnostic classification and its analysis For the sought P/D area ratios of all investigated 240 healthy eye images and 240 glaucomatous eye images, the ð2Þ ð240Þ corresponding real neural net outputs yð1Þ h ; yh ; . . . ; y h ð1Þ ð2Þ ð240Þ and yg ; yg ; . . . ; yg were computed. The obtained histograms of the statistical distribution of neural net outputs yh in the case of healthy eye images and of the statistical distribution of neural net outputs yg in the case of the glaucomatous eye images are shown in Figs. 6 and 7. The eyes, investigated by the described neural net method, were sorted out into the three classes: lowglaucoma-risk class, glaucoma-suspected class and highglaucoma risk class. The low-glaucoma-risk class includes eyes with the corresponding neural net outputs below t1, the glaucoma-suspected class involves eyes with the corresponding neural net outputs within the interval ht1 ; t2 i and the high-glaucoma-risk class includes eyes with the corresponding neural net outputs above t2. At the same time, the parametric values t1 and t2, 0ot1 ot2 o1; are the suitable distinguishing limits. The chosen values t1 ¼ 0:13; t2 ¼ 0:63 are considered in this chapter. Thus the probability Ph/high, that the healthy eye is inconveniently classified to the high-glaucoma-risk class, is Ph=high ¼ 1 Dh ðt2 Þ ¼ 0:05
(10)
and the probability Pg/low, that the glaucomatous eye is inconveniently classified to the low-glaucoma-risk
Fig. 6. Obtained histogram of the statistical distributions of the neural net outputs yh in the case of healthy eye images.
ARTICLE IN PRESS F. Pluha´cˇek, J. Pospı´ sˇ il / Optik 116 (2005) 133–143
139
which expressions follow from the common errorpropagation rule of the explanation in [17]. The concrete values sh=high 0:01; sg=low 0:02; sh=suspected 0:04; sg=suspected 0:02 were numerically obtained by computation for the considered distinguishing values t1 ¼ 0:13; t2 ¼ 0:63; marked in Figs. 6 and 7.
6. Probability density curves method for diagnostic evaluation of the total P/DT area ratios 6.1. Probability densities estimation of the P/DT statistical distributions
Fig. 7. Obtained histogram of the statistical distributions of the neural net outputs yg in the case of glaucomatous eye images.
The probability densities ph and pg of the experimentally obtained statistical distributions of the x ¼ P=DT values of the above-mentioned 240 healthy and 240 glaucomatous OD images can be approximated by the experimentally ascertained relations
class, is
ph ðP=DT Þ ¼ N h exp½bh1 ðP=DT bh2 Þ3 ; P=DT 2 h0; 1i; (18)
Pg=low ¼ Dg ðt1 Þ ¼ 0:05:
(11)
The above-introduced functions Dh(t) and Dg(t) are experimentally obtained statistical distribution functions of the neural net outputs yh and yg. On the other hand, the probability Ph/suspected, respectively, Pg/suspected, for the classification of healthy, respectively, glaucomatous, eye to the glaucoma-suspected class can be computed by the relations Ph=suspected ¼ Dh ðt2 Þ Dh ðt1 Þ ¼ 0:45;
(12)
pg ðP=DT Þ ¼ N g exp½bg1 ðP=DT bg2 Þ2 ; P=DT 2 h0; 1i; (19) where Z
1
Nh ¼ 1
exp½bh1 ðx bh2 Þ3 dx;
(20)
exp½bg1 ðx bg2 Þ2 dx
(21)
0
Z
1
Ng ¼ 1 0
Pg=suspected ¼ Dg ðt2 Þ Dg ðt1 Þ ¼ 0:34:
(13)
With respect to the already mentioned non-zero mean standard deviation sn of the neural net outputs, the introduced probabilities Ph/high, Pg/low, Ph/suspected and Pg/suspected were determined with the relevant standard deviations @Ph=high (14) sn ; sh=high @t2 @Pg=low sg=low (15) sn ; @t1 " sh=suspected sn
@Ph=suspected @t1
2
#1=2 @Ph=suspected 2 þ ; @t2
are the real normalization constants and bh1 ; bh2 ; bg1 ; bg2 (shortly denoted by bh1;2 and bg1;2 ) represent the suitable real parameters. The functions (18) and (19) were found by histograms of P/DT values of healthy and glaucomatous eye images using the least-squares mathematical method described in [17,18]. For each histogram, the 10 class intervals were utilized. All computations were performed numerically with a sufficient accuracy. The obtained estimations of parameters bh1;2 ; bg1;2 and their standard deviations sh1 ; sh2 ; Table 1. Using least-squares method obtained estimations of parameters bh1;2 ; bg1;2 and their standard deviations sh1;2 ; sg1;2 Parameter
Estimated value
Standard deviation
Estimated value
bh1 bh2 bg1 bg2
7.0 102 11.4 102 5.2 101 24.3 102
sh1 sh2 sg1 sg2
0.6 102 0.4 102 0.4 101 0.4 102
(16) sg=suspected
"
#1=2 @Pg=suspected 2 @Pg=suspected 2 sn þ ; @t1 @t2 (17)
ARTICLE IN PRESS 140
F. Pluha´cˇek, J. Pospı´ sˇ il / Optik 116 (2005) 133–143
Fig. 8. Resultant histogram of the statistical distribution of the healthy eye P/DT images with the probability density curve ph(P/DT) of expression (18).
Fig. 9. Resultant histogram of the statistical distribution of the glaucomatous eye P/DT images with the probability density pg(P/DT) of expression (19).
sg1 ; sg2 (shortly expressed by sh1;2 and sg1;2 ) are introduced in Table 1. Both these probability densities (18) and (19) were also tested by the chi-squared test of fit [17,18] with the result that they did not reject for the significant levels below 22% in the case of ph and below 42% in the case of pg. The obtained resultant histograms of the healthy and glaucomatous eye images P/DT with the corresponding curves of the probability densities ph, pg are presented in Figs. 8 and 9.
6.2. Proposal of the probability density curves for the eye diagnostic classification and its analysis For the eye diagnostic classification, based on the described probability density curves method, the following three glaucoma risk classes of eyes were considered in the modified analogy with the previous neural net method: low-glaucoma-risk class, glaucomasuspected class and high-glaucoma-risk class. The
Fig. 10. Graphical dependences of the probability densities ph and pg on the P/DT area ratios under validity of the relations (18) and (19) for the in Table 1 introduced values of bh1;2 and bg1;2 :
low-glaucoma-risk class includes the eyes with the P/DT below t1, the glaucoma-suspected class involves the eyes with the P/DT from the interval ht1 ; t2 i and the high-glaucoma-risk class includes the eyes with the P/DT above t2. The distinguishing parametric limits t1 and t2, 0ot1 ot2 o1 are the suitably chosen values of P/DT, see Fig. 10. Under these conditions, the healthy eye is inconveniently classified to the high-glaucoma-risk class with the probability Z 1 Ph=high ¼ ph ðxÞ dx; (22) t2
while the glaucomatous eye is inconveniently classified to the low-glaucoma-risk class with the probability Z t1 Pg=low ¼ pg ðxÞ dx: (23) 0
With respect to the nonzero values of the standard deviations sT ; sh1;2 and sg1;2 ; the introduced probabilities Ph/high and Pg/low were determined with the relevant standard deviations sh=high ; sg=low of approximate definitions @Ph=high sT ; (24) sh=high @t2 @Pg=low sg=low sT ; (25) @t1 which follow from the common error-propagation rule. The performed numerical analysis showed that the influence of sh1;2 and sg1;2 on sh=high ; sg=low is negligible and therefore need not be considered. For the chosen values t1 ¼ 0:09; t2 ¼ 0:22; the equalities Ph=high ¼ 0:05; sh=high 0:02; Pg=low ¼ 0:05; sg=low 0:01 are valid.
ARTICLE IN PRESS F. Pluha´cˇek, J. Pospı´ sˇ il / Optik 116 (2005) 133–143
141
The remaining probability Ph/suspected for the classification of a healthy eye in the glaucoma-suspected class and the remaining probability Pg/suspected for the classification of a glaucomatous eye in the glaucomasuspected class now fulfil the relations Z t2 Ph=suspected ¼ ph ðxÞ dx; (26) t1
Z
t2
pg ðxÞ dx:
Pg=suspected ¼
(27)
t1
The corresponding relevant standard deviations sh=suspected ; sg=suspected for Ph/suspected and Pg/suspected can be computed using again the above-mentioned common error-propagation rule, which gives the equations "
#1=2 @Ph=suspected 2 @Ph=suspected 2 sh=suspected sT þ ; @t1 @t2 (28)
sg=suspected
"
#1=2 @Pg=suspected 2 @Pg=suspected 2 sT þ : @t1 @t2 (29)
In our case for the considered equalities t1 ¼ 0:09 and t2 ¼ 0:22; the values Ph=suspected ¼ 0:45; sh=suspected 0:04; Pg=suspected ¼ 0:34 and sg=suspected 0:03 were obtained.
6.3. Glaucoma prevalence probability in view of the total P/DT area ratios and the subject ages For the total probability P of the occurrence of glaucomatous eyes in a population group of eyes with the given P/DT value and its standard deviation sT (in short called the glaucoma prevalence probability), the corresponding formula P¼
P p pg Pp pg þ ð1 Pp Þph
(30)
was derived on the basis of the conditional probability theory [18]. In our case, the used probability densities ph and pg are given by the relations (18) and (19). The introduced parameter Pp with the standard deviation sPp means the probability of glaucoma occurrence in the whole population group (without reference to the P/DT values) of the given age range. In general, Pp is markedly influenced by the race of the investigated subjects and differs for different age ranges, see [1,21]. The relevant resultant standard deviation sP for the total probability P was estimated according to the common errorpropagation rule by the suitable relation
Fig. 11. Graphs of the quantities (30) and (31) versus the P/DT area ratios constructed for the examined subject age ranges below 65 years and above 75 years in the case of the parametric values of bh1;2 ; bg1;2 and sh1;2 ; sg1;2 introduced in Table 1.
" sP
@P @ðP=DT Þ
2
s2T
@P þ @Pp
2
#1=2 s2Pp
;
(31)
where the influence of sh1;2 and sg1;2 was acceptably neglected. Relations (30) and (31) are plotted as functions of P/DT in Fig. 11 for two different age ranges of the white population (below 65 years and above 75 years) and for in Table 1 introduced estimated values of bh1;2 ; bg1;2 ; sh1;2 and sg1;2 : In the case of age below 65 years, the values Pp ¼ 0:01; sPp ¼ 0:005 were obtained and for age above 75 years, the equalities Pp ¼ 0:03; sPp ¼ 0:01 are valid, see [21]. It is evident from the graphs in Fig. 11 that the glaucoma prevalence probability P increases with the increased P/DT values in both considered cases. Moreover, the resultant standard deviations sP show the maximum and their values are inconveniently large, mainly for the interval of the maximal slope of P.
7. Comparison of the presented methods and the obtained results The described quantile curves method, based on the utility of the quantile standpoint for diagnostic evaluation of the healthy and glaucomatous human eyes, enables advantageously the convenient direct objective determination of the corresponding quantile curves and then also of the number of the representative OD sectors with undesirably increased values of the P/D area ratios. The disadvantage of this method is that the final glaucoma diagnosis must be often accompanied with a subjective decision without any possibility to compute directly the resultant diagnostic probability. On the other hand, the also presented neural net and probability density curves methods enable
ARTICLE IN PRESS F. Pluha´cˇek, J. Pospı´ sˇ il / Optik 116 (2005) 133–143
142
advantageously an objective final automatic classification of the investigated eyes in three exactly defined glaucoma risk classes. Because probability values Ph/high or Pg/low for classification of healthy or glaucomatous eye in the high-glaucoma or low-glaucoma risk class are conveniently small enough for both later-mentioned glaucoma diagnostic methods (Ph=high ¼ Pg=low ¼ 0:05), the eyes from the high-glaucoma risk class, respectively, low-glaucoma risk class, can be reliably diagnosed as really glaucomatous, or healthy. In such cases, the probability of incorrect diagnosis is namely also equal to Ph/high, or Pg/low, i.e. it is small. The remaining probabilities Ph/suspected and Pg/suspected for classification of healthy or glaucomatous eye in the glaucomasuspected class are both larger (Ph=suspected ¼ 0:45 and Pg=suspected ¼ 0:34). Thus a reliable diagnosis of the eye must be carried out differently, by means of a characteristic different quantity. According to Machala and Pospı´ sˇ il [22], the suitable reliability factor of possible definition r¼
j¯th ¯tg j ½ðs2h
þ s2g Þ=21=2
(32)
can be applied. The corresponding characteristic parameters ¯th ; ¯tg and sh ; sg for healthy and glaucomatous eyes are the mean values and standard deviations relating to the obtained statistical distribution functions Dh(t) and Dg(t) for the outputs t ¼ yh and yg in the case of the neural net method, or the mean values and standard deviations of the obtained distributions of the probability densities ph(t) and pg(t) for the values t ¼ P=DT in the case of the probability density curves method. From the definition (32) it is evident that the value of r; representing in fact the normalized distance between ¯th and ¯tg ; increases with increased difference j¯th ¯tg j and decreases with increased standard deviations sh ; sg : A higher value of r corresponds to a higher diagnostic reliability of the used method and vice versa. Values r ¼ 1:82 and 1.72 were obtained for the neural net method and the probability density curves method. Because these values differ only little for both considered methods, they have an approximately identical reliability. In addition, the probability density curves method enables the computation of the glaucoma prevalence probability for the given total P/DT area ratios and the subject age ranges. With respect to these facts and to the complicated structure of the utilized neural net, the probability density curves method is mostly more useful in practical applications.
8. Conclusion The contribution of the present article consists in the proposal, verification and comparison of three new
computer image analysis statistical methods for experimental diagnostic evaluation of colour glaucomatous changes within the optic disc (OD) of a human eye retina. The search for such changes by means of a suitable computer-connected digital fundus camera is based on the preliminary determination of the total representative OD areas and their relevant relative pallor area sizes (P/D), which are represented by 36dimensional characteristic vectors of components xn ¼ P=D1 ; P=D2 ; . . . ; P=D36 or, respectively, by five-dimensional characteristic vectors of components xk ¼ P=DT ; P/DI, P/DII, P/DIII, P/DVI. The suitable established components relate to the elementary P/D area ratios for 36 equally spaced OD sectors or, respectively, to the combination of the total OD area (P/DT) and the partial quadrants numbered I, II, III and IV in the total OD area for the sufficiently large representative sets of the normal healthy and the glaucomatous eyes. As the first is presented the so-called the quatnile curves method, based on the knowledge of the vector components P/Dn, n ¼ 1; 2; . . . ; 36; and their experimentally obtained distribution functions Dhn(P/Dn) and Dgn(P/Dn) for the sets of 240 normal healthy and 240 glaucomatous eyes and for the special parametric values ðxÞ qðxÞ hn and qgn of P/Dn relating to the suitable chosen normal and glaucomatous quantile x: The resultant main dependences of these well-arranged values in the descending order on their sequence in this ordering represent the x% normal and glaucomatous quantile curves for the successive and effective evaluation purposes also in comparison with the suitably introduced one-eye normal and glaucomatous curves. The second investigated method, specified as the neural net method, utilizes the modified artificial computer-simulated neural net point of view with the corresponding back-propagation learning algorithm for diagnostic evaluation of the 240 normal healthy and 240 glaucomatous eyes on the basis of their fivedimensional characteristic input vectors of type x ¼ ðP=DT ; P=DI ; P=DII ; P=DIII ; P=DIV Þ: For the searched P/D area ratios of all investigated healthy and glaucomatous eye images, the corresponding resultant neural ð2Þ ð240Þ ð2Þ ð240Þ net outputs yð1Þ and yð1Þ g ; yg ; . . . ; yg h ; yh ; . . . ; yh were computed and their histograms were constructed. On the basis of the resultant histograms, the investigated eyes were classified into three glaucoma risk classes and diagnosed by means of the significant error probabilities and standard deviations. The last considered method, i.e. the probability density curves method, relates to the obtained statistical distributions of the P/DT values for the also considered 240 healthy and 240 glaucomatous OD images. It is based on the exploitation of the experimentally ascertained relations for the probability density curves ph(P/DT) and pg(P/DT) relating to the mentioned sets of the investigated normal healthy and glaucomatous
ARTICLE IN PRESS F. Pluha´cˇek, J. Pospı´ sˇ il / Optik 116 (2005) 133–143
eyes and to a suitably chosen parameters. These curves were found by the corresponding histograms of the P/ DT values using the least-squares mathematical method. On the basis of these curves, the investigated eyes were classified into three glaucoma risk classes, on analogy with the previous neural net method, but by means of the significant error probabilities and standard deviations of different concrete definitions. Moreover, also the total glaucoma prevalence probability and the corresponding resultant standard deviation were estimated for different age ranges of the investigated subjects. Both the last two automatic and objective statistical diagnostic methods can replace effectively a subjective evaluation of glaucomatous changes within the representative OD. The performed analysis shows that the reliabilities of these methods are approximately identical. In comparison with other classical non-statistical glaucoma diagnostic methods, based directly on the numerical representation of the OD glaucomatous changes (e.g. the stereoscopic images computer evaluation method or the laser scanning tomography method), the designed statistical methods require only the computer-connected digital fundus camera arrangement without other basic devices. They enable also a larger diversity of the current ophthalmologic equipment without any new expensive investment.
Acknowledgement All the investigated optic disc images were kindly supplied by the physician MUDr. Toma´sˇ Kubeˇna. The authors would like to express their gratitude to him.
References [1] B. Shields, Textbook of Glaucoma, Williams & Wilkins, Baltimore, 1992. [2] H. Kraus, J. Boguszakova´, P. Diblı´ k, D. Dotrˇ elova´, M. Filipec, J. Hycl, I. Karel, I. Kocour, H. Kraus, J. Krˇ epelkova´, P. Kuchynka, J. Otradovec, J. Palecˇkova´, J. Peregrin, J. R˘ehu˚rˇ ek, J. Sveˇra´k, A. Sˇedivy`, Compendium of the Eye Medicine, Grada, Praha, 1997 (in Czech). [3] A.P. Nesterov, Primary Glaucoma, Avicenum, Praha, 1991 (in Czech). [4] S. R˘eha´k, R. Knobloch, O. Riebel, F. Vrabec, Medicine of the Eye, Avicenum, Praha, 1980 (in Czech). [5] A. Nover, The Ocular Fundus. Methods of Examination and Typical Findings, Schattauer, Stuttgart, 1987.
143
[6] H. Sauter, W. Straub, H. Rossmann, Atlas of the Ocular Fundus. Photographs of Typical Changes in Ocular and Systemic Disease, Urban & Schwarzenberg, Mu¨nchen, 1997. [7] T.N. Cornsweet, S. Hersh, J.C. Humphries, R.J. Beesmer, D.W. Cornsweet, Quantification of the shape and color of the optic nerve head, In: G.M. Breinin, I.M. Siegel (Eds.), Advances in Diagnostic Visual Optics, Springer, New York, 1983, pp. 141–149. [8] R. Varma, G.L. Spaeth, The par is 2000: a new system for retinal digital image analysis, Ophthal. Surg. 19 (1988) 183–192. [9] T.N. Cornsweet, S. Hersh, Ocular-fundus analyzer. US Patent No. 4715703, US Government, Washington, 1987. [10] F. Pluha´cˇek, J. Pospı´ sˇ il, Present methods of the computer-image analysis of a human retina with regard to diagnostics of glaucoma, Fine Mech. Opt. 46 (2001) 326–334 (in Czech). [11] T. Shioiri, Ophthalmologic image processor, US Patent No. 6104828, US Government, Washington, 2000. [12] F.S. Mikelberg, K. Wijsman, M. Schulzer, Reproducibility of topographic parameters obtained with the heidelberg retina tomograph, J. Glaucoma 2 (1993) 101–103. [13] G. Zinser, U. Harbarth, H. Schro¨der, Formation and analysis of three-dimensional data with the laser tomographic scanner (LST), In: J.E. Nasemann, R.O.W. Burk (Eds.), Scanning Laser Ophthalmoscopy and Tomography, Quintessenz Verlag, Mu¨nchen, 1990, pp. 243–252. [14] K.U. Bartz-Schmidt, A. Sengersdorf, P. Esser, P. Walter, R.D. Hilgers, G.K. Krieglstein, The cumulative normalised rim/disc area ratio curve, Greafe’s Arch. Clin. Exp. Ophthalmol. 234 (1995) 227–231. [15] J. Caprioli, Reproducibility of optic disc measurements with computized analysis of stereoscopic video images, Arch. Ophthalmol. 104 (1986) 1035–1039. [16] T. Kubeˇna, P. Cˇernosˇ ek, J. Mayer, Technique of digital planimetry of optic disc, Czech Slovak Ophthalmol. 56 (2000) 170–175 (in Czech). [17] L. Kuba´cˇek, L. Kuba´cˇkova´, Statistics and Metrology, Palacky` University Press, Olomouc, 2000 (in Czech). [18] J. Sˇkra´sˇ ek, Z. Tichy`, Principles of Applied Mathematics III, SNTL, Praha, 1990 (in Czech). [19] B. Mu¨ller, J. Reinhardt, M.T. Strickland, Neural Networks: An Introduction, Springer, Berlin, 1995. [20] J.G. Taylor, Neural Networks and their Applications, Wiley, Chichester, 1996. [21] B.E. Klein, R. Klein, W.E. Sponsel, T. Franke, L.B. Cantor, J. Martone, M.J. Menage, Prevalence of glaucoma: the beaver dam eye study, Ophthalmol. 99 (1992) 1499–1504. [22] L. Machala, J. Pospı´ sˇ il, Proposal and verification of two methods for evaluation of the human iris video-camera images, Optik 112 (2001) 335–340.
Optik
Optics
Optik 116 (2005) 133–143 www.elsevier.de/ijleo
Proposal, verification and comparison of three computer image analysis methods for detection and evaluation of colour glaucomatous changes within the optic disc of a human eye retina Frantisˇ ek Pluha´cˇeka,, Jaroslav Pospı´ sˇ ilb a
Department of Optics of Palacky` University, 17. Listopadu 50, 77200 Olomouc, Czech Republic Department of Experimental Physics and Joint Laboratory of Optics of Palacky` University and Institute of Physics of Academy of Sciences, 17. Listopadu 50a, 77200 Olomouc, Czech Republic
b
Received 20 August 2004; accepted 28 November 2004
Abstract The typical symptom of the human eye glaucoma is a rise and a progression of the bright area (named pallor area) within the retina blind spot. The image analysis manner proposed by the authors detects and suitably numerically describes the relative size of the representative pallor area in the colour digital image of the retina obtained by a suitable fundus camera connected with the computer. Three new different computer image analysis statistical methods for experimental diagnostic evaluation of the obtained characteristic data are proposed in this article: the quantile curves method, the neural net method and the probability density curves method. The quantile curves method is based on the graphical comparison of a relative representative pallor area size with its determined normal value. The neural net and probability density curves methods can automatically and objectively classify the investigated eyes in exactly defined glaucoma risk classes and diagnosed glaucoma with the rated probabilities of incorrect diagnosis determination. All mentioned methods are verified and mutually compared by their application to the large statistical sets of human retina images of various healthy and glaucomatous subjects. r 2005 Elsevier GmbH. All rights reserved. Keywords: Glaucoma; Optic disc; Image analysis; Statistical evaluation; Neural net
1. Introduction Glaucoma is the common name of a group of single human eye diseases characterized by an increased intraocular pressure (the eyewater pressure in the eyeball), typical damage to the photosensitive cell nerve fibres of the retina, colour and shape changes of the blind spot of the retina (named optic disc, OD) and in Corresponding author. Tel.: +420 585634310.
E-mail addresses: [email protected] (F. Pluha´cˇek), [email protected] (J. Pospı´ sˇ il). 0030-4026/$ - see front matter r 2005 Elsevier GmbH. All rights reserved. doi:10.1016/j.ijleo.2004.11.009
the eye visual field [1–6]. The prevalence of glaucoma is approximately in 1–2% of adults [1]. The patient with this disease feels no subjective problems for a long time, but there is a slow decrease of the visual functions. The final stage of non-medicate glaucoma is practical blindness. Thus, the early diagnostics of glaucoma is very important. For the progression of this disease, the size of an iridocorneal angle can be representative, see publications [1–4]. The considered iridocorneal angle is the anterior chamber space between the iris base and the inner cornea margin. Its size influences the intraocular pressure and the glaucoma risk. In view of the discussed
ARTICLE IN PRESS 134
F. Pluha´cˇek, J. Pospı´ sˇ il / Optik 116 (2005) 133–143
angle size, two glaucoma types are commonly distinguished: the open-angle glaucoma (if the angle is of normal size) and the angle-closure glaucoma (if the angle is close). Only the routine open-angle glaucoma is considered in this article. One of the first symptoms of glaucoma is the decay of the retina photosensitive cell nerve fibres and a subsequent progress of the depression or cup within the OD. These glaucomatous processes cause a rise and enlargement of the bright area (named pallor area) within the OD cup bottom, which represents the OD area with the absence of the nerve fibres [1]. In the case of a healthy eye, the pallor area is missing or is only small. Because the subjective determination and evaluation of the pallor area or other symptoms have small reproducibility, the objective diagnostic methods are required. The size of the cup can be objectively measured and evaluated at present by the method of the computer evaluation of stereoimages [7–11] and by the laser scanning tomography method [12–15]. It is also possible to scan the human eye retina using a suitable colour fundus camera connected to the computer by a convenient graphic card, as was carried out by the authors of this article. A colour digital image obtained in this way can be advantageously archived in an appropriate data medium and can be used for detection of the colour glaucomatous changes using a suitable computer (digital) image analysis. For diagnostic purposes, the relative size of the so-called representative pallor area in the image must be considered because of its dependence on the considered so-called representative OD size. It is the ratio of the representative pallor area size and the representative OD area size in the image (denoted as P/D area ratio or only P/D) which is chosen as the relative pallor area size measure in the case of the utilized image analysis manner presented in this article. The representative OD area is conventionally divided into several suitable parametric sectors and the P/D area ratio is computed for each sector. For diagnostic evaluation of the P/D area ratios, three new statistical methods of computer image analysis are proposed, verified and compared in this paper. In accordance with the utilized principle and form of the obtained main results, the following names for these methods were accepted: the quantile curves method, the neural net method, and the probability density curves method. All computer calculations were realized by the mathematical programs MATLAB and MATHCAD.
2. Used equipment and the measuring conditions The colour computer digital OD images of the size 576 768 pixels, obtained by the computer-connected
digital fundus camera DFK 98 [16], were used. The mentioned camera contains the colour CCD videocamera connected with the adapted photographic system Retinofot Carl Zeiss Jena 201, see [16]. During the photographic scanning of the OD, the described equipment was located so that the rows of the acquired image matrix correspond to the horizontal direction. The actual size of every OD area was approximately 5.8 mm 4.4 mm2. The brightness values of each of the corresponding three image colour channels (i.e. red, green and blue channel) were quantized to 256 discrete levels and represented by integers from 0 to 255. Value 0 represents the minimum brightness and value 255 the maximum brightness of each colour digital image. All the following considerations and results relate to these measuring conditions. The maximum number of 240 OD images of 120 normal healthy subjects and of 240 OD images of 120 subjects with the open-angle glaucoma were utilized in our research. According to an eventual medical requirement, additional medical diagnoses were proposed by the cooperate glaucomaspecialized ophthalmologist on the basis of independent measurement of the intraocular tension, visual field, enlargement of the cup using the laser scanning tomograph and of a subjective observation of the thickness of the retina nerve fibre layer. The examined normal healthy subjects had the age range from 14 to 76 years with the average age of 46 years and with the standard deviation of 13 years. The glaucomatous subjects under examination had the age range of 10–82 years with the average of 51 years and the standard deviation of 16 years.
3. Performed preliminary determination of the representative OD area and the relative pallor area size by the utilized computer image analysis manner 3.1. Brightness of the representative OD area and definition of the representative pallor area The representative pallor area is situated within the representative OD area in the retina image obtained by the above-mentioned fundus camera. The representative OD area was manually demarcated by a convenient little excentric representative boundary ellipse, which connects the five auxiliary points selected by the operator’s experience in the used colour digitized image of the OD. In order to meet the requirement of objective orientation of the image, a polar coordinate system was chosen with its pole 0 in the centre of the representative boundary ellipse and with the polar line oriented horizontally according to the reverse direction of the human nose, see
ARTICLE IN PRESS F. Pluha´cˇek, J. Pospı´ sˇ il / Optik 116 (2005) 133–143
Fig. 1. The angles were considered clockwise for the right eye and anticlockwise for the left eye. It was ascertained experimentally that the OD image brightness G is most expressive in the green colour channel and is displayed with sufficient high optical contrast. Hence, only this channel was used for the needed normalized brightness determination. For such normalization the two suitable reference pixel groups A and B were considered in accordance with Fig. 1. It is evident that the image point (u,v) located in the intersection of the uth row and the vth column of the pixel array lies within these pixel groups of empirically chosen width a/5 when hold the inequalities qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a rðfÞ p ðu u0 Þ2 þ ðv v0 Þ2 prðfÞ; (1) 5 where 15 pfp 15 in the case of the group A or 165 pfp195 in the case of the group B. The point ðu0 ; v0 Þ represents the position of the pole 0 of the abovedefined polar coordinate system, a is the reference major semi-axis size of the representative boundary ellipse, f represents the oriented angle between the position vector of a pixel (u,v) and the polar line, rðfÞ is the mutual distance between the pole 0 and the boundary ellipse point of the angular coordinate f: Because the brightness values of the OD image areas corresponding to the pixel groups A, B are not practically influenced by the glaucoma progression, they can be used to determine the reference brightness values GA and GB for the already mentioned normalization of G. These values are defined as the green image channel brightness values of the largest occurrence in pixels of the groups A and B. Moreover, the location of each pixel group can be represented by the average values u¯ A ; v¯ A and u¯ B ; v¯ B : The also required normalized brightness Gn(u,v) in each point (u,v) of the representative OD image is proposed by the useful relation G n ðu; vÞ ¼
2Gðu; vÞ : GA þ GB
(2)
135
Fig. 2. Example of the healthy (left) and glaucomatous (right) OD image with the corresponding negligible small and expressive large bright pallor area. The boundary ellipses define the representative OD areas.
The above-mentioned values GA and GB differ in general, because the brightness of the OD periphery often inconveniently changes in the polar line direction. For the effective compensation of this effect, we designed and experimentally verified the acceptable brightness transformation G 0n ðu; vÞ ¼ G n ðu; vÞ G An
G Bn G A n ðv v¯ A Þ v¯ B v¯ A
(3)
of the normalized brightness Gn(u,v) in the pixel (u,v) to the new value G0n ðu; vÞ: The introduced values G An and G Bn are the values GA and GB, which are normalized by value (2). The corresponding representative pallor area of the OD image is then defined as the image area within the representative boundary ellipse, where the normalized transformed value G 0n is greater than or equal to the selected suitable threshold value p. The empirically established equality p ¼ 0:2 is accepted in the following text. Fig. 2 shows an example of healthy (left) and glaucomatous (right) OD image with the corresponding negligible small and expressive large bright pallor area. The boundary ellipses define the representative OD areas.
3.2. Relative pallor area size characteristic vectors and the standard deviations of their measurement
Fig. 1. Example of an optic disc representative boundary ellipse, the considered polar coordinate system and the chosen pixel groups A and B.
Let us assume the knowledge of the representative OD and pallor areas which were determined in agreement with the previous text. The relevant relative pallor area size P/D, mentioned in Chapter 1 of this article and related to the total representative OD area, can be conventionally numerically represented by the 36-dimensional characteristic vector (P/D1, P/D2,y, P/D36) or by the five-dimensional characteristic vector (P/DT, P/DI, P/DII, P/DIII, P/DIV) in accordance with [14,15]. The introduced single components P/D1, P/D2,y,P/D36 represent the elementary P/D area ratios for 36 equally spaced OD sectors with the common vertex in the pole 0 of the used polar coordinate system
ARTICLE IN PRESS 136
F. Pluha´cˇek, J. Pospı´ sˇ il / Optik 116 (2005) 133–143
and with the vertex angle ranges h0 ; 10 Þ; h10 ; 20 Þ; . . . ; h350 ; 360 Þ: On the other hand, the single components of the five-dimensional characteristic vector (P/DT, P/DI, P/DII, P/DIII, P/DIV) denote successively the P/D area ratios of the total representative OD area and of its four quadrants numbered I, II, III, IV and contained the vertex angle ranges h45 ; 45 Þ; h45 ; 135 Þ; h135 ; 225 Þ; and h225 ; 315 Þ: The performed experiments showed that the used subjective manual demarcation of the representative OD elliptical border for a given OD image is markedly influenced by a significant variability, which can be represented by the single standard deviations of the above defined and measured characteristic components. For the performed estimation of such measurement single deviations, the demarcation of the every representative OD elliptical border and also the subsequent computation of all the corresponding partial P/D area ratios were always repeated ten times for randomly selected 40 healthy and 40 glaucomatous OD images (i.e. 800 demarcations were performed). Then for each selected image, the obtained partial P/D area ratios were statistically processed in the way that the associated single standard deviation ss ; s 2 fT; I; II; III; IV; 1; . . . ; 36g; of the repeated measurement of a partial P/Ds area ratio is represented by the mean value of all corresponding partial P/Ds standard deviations. Furthermore, for practical usage, the single values s1 ; s2 ; . . . ; s36 can be often suitably replaced by their average of the possible expression s1;2;...;36 : Our computations yield the values sT ¼ 0:008; sI ¼ 0:02; sII ¼ 0:01; sIII ¼ 0:007; sIV ¼ 0:01; s1;2;...;36 ¼ 0:02: The results of investigated dependences of the P/DT discrete values on the age of the examined 240 healthy and 240 glaucomatous subjects were suitably approximated by the corresponding regression lines using the mathematical theory of regression analysis [17,18]. The appropriate regression lines of equations P=DT ¼ 5:6 105 age þ 0:1; P=DT ¼ 1:7 103 age þ 0:3 with the correlation coefficients 0.01, 0.3 and with the standard deviation 0.3 for both regression line slopes were obtained by computation. The used variable quantity under indication age means the subject age in years. The obtained small values of the regression line slopes and the relative high values of their standard deviations indicate that the age dependences are insignificant for the considered subject sets.
4. Quantile curves method for diagnostic evaluation of the 36-dimensional P/D area ratios The above-defined 36-dimensional characteristic vector (P/D1, P/D2,y,P/D36) of the P/D area ratio, which
is fundamental for this method, was computed for each of the 240 healthy and 240 glaucomatous OD images mentioned in Chapter 2. The obtained statistical probability distributions of the values of each single component P/Dn, n ¼ 1; 2; . . . ; 36; can be advantageously mathematically compared by the associated experimental distribution functions Dhn(P/Dn) in the case of healthy OD images and Dgn(P/Dn) in the case of glaucomatous OD images. According to [14], the special ð90Þ ð50Þ ð90Þ parametric values qð50Þ of P/Dn, hn ; qhn and qgn ; qgn fulfilling the general relations Dhn ðqðxÞ hn Þ ¼ x=100;
x 2 f50; 90g;
(4)
Dgn ðqðxÞ gn Þ ¼ x=100;
x 2 f50; 90g
(5)
and named as x% normal quantiles and x% glaucomatous quantiles, were utilized. Thus, the probability of occurrence of the P/Dn values in normal healthy eyes, or glaucomatous eyes, up to the extreme value ðxÞ qðxÞ hn ; respectively, qgn is x%: For our comparison, the ðxÞ ðxÞ ðxÞ ðxÞ ðxÞ considered values qðxÞ h1 ; qh2 ; . . . ; qh36 and qg1 ; qg2 ; . . . ; qg36 of P/Dn were at first well-arranged in the descending order and then plotted in dependence on their sequence in this ordering (see Fig. 3). The obtained graphs represent the x% normal quantile curves (solid lines) and the glaucomatous quantile curves (dashed lines) for the considered sets of eyes. Moreover, the 36-dimensional characteristic vector can be also applied only to one arbitrary eye. Such eye is graphically described by a one-eye curve which represents the dependence of its single P/Dn values sorted in the descending order on their sequence in this ordering. Examples of the
Fig. 3. The obtained 50%- and 90%- normal quantile curves (solid lines) in comparison with the 50%- and 90%glaucomatous quantile curves (dashed lines) and with oneeye curves of one randomly selected healthy (normal) eye (dotted line) and glaucomatous eye (dot–dashed line). These curves relate to the described quantile curves method.
ARTICLE IN PRESS F. Pluha´cˇek, J. Pospı´ sˇ il / Optik 116 (2005) 133–143
137
obtained one-eye curves for one randomly selected healthy (normal) eye (dotted line) and glaucomatous eye (dot–dashed line) are also shown in Fig. 3. It is evident that the values of P/D area ratios, presented in the graphical form for the glaucomatous eye, are markedly influenced by the glaucoma. Generally, if the one-eye curve of an investigated eye is above the 90% normal quantile curve, the high risk of glaucoma occurrence in this eye can be assumed. Such a graphical comparison can be a useful glaucoma diagnostic tool.
5. Neural net method for diagnostic evaluation of the five-dimensional P/D area ratios 5.1. Description of the structure of the used artificial computer-simulated neural net Generally, the artificial neural net represents the mathematical transformation system with a N-dimensional vector x ¼ ðx1 ; . . . ; xN Þ of input components and a M-dimensional vector y ¼ ðy1 ; . . . ; yM Þ of output components, see publications [19,20]. The basic element of every neural net is a formal neuron, whose inputs contain the information from other neurons or from the outer environment. Each neuron has only one output, transmitting the resulting information from this neuron to further elements (neurons) or to the output of the net. All the elements are arranged in L layers and each neuron in the lower layer is connected with each neuron in the layer above it. The neurons in the same layer are not mutually connected. Inputs of the neurons of the lowest layer (denoted as the first or input layer) are real inputs of the net and outputs of the neurons of the highest layer (denoted as the Lth or output layer) represent the real outputs of the net. According to Refs. [19,20], the mathematical formal model of the kth neuron in the ith layer can be represented by the equation " !#1 N i1 X xk;i ¼ 1 þ exp wi1;j;k xj;i1 yk;i ; (6) j¼1
where xk,i is the output of the kth neuron in the ith layer, xk;0 ¼ xk is the kth input of the whole neural net, xk;L ¼ yk is the kth output of the whole neural net, wi,j,k represents the weight of connection of the jth neuron in the ith layer with the kth neuron in the (i+1)th layer and w0,j,k is the weight of connection of the jth neural net input with the kth neuron in the first layer ði ¼ 1Þ: Furthermore, Ni denotes the number of neurons in the ith layer, N 0 ¼ N is the number of the neural net inputs (called the dimension of the input), N L ¼ M is the number of the neurons in the output layer (called the dimension of the output) and yk;i represents the
Fig. 4. Scheme of the considered neural net structure for diagnostic evaluation of the five-dimensional P/D area ratios.
threshold of the kth neuron in the ith layer. It comes true that the output signal of each neuron in the ith layer represents the input of all neurons of the (i+1)th layer at the same time. The neural net with the parameters L ¼ 2; N ¼ N 0 ¼ 5; N 1 ¼ 3; M ¼ N 2 ¼ 1 and with the input vector x ¼ ðP=DT ; P=DI ; P=DII ; P=DIII; P=DIV Þ is considered in the following text. Its structural scheme is shown in the Fig. 4.
5.2. Back-propagation learning algorithm The values of the already mentioned thresholds and weights of the considered neurons in the neural net were determined by the so-called back-propagation learning algorithm, see [19,20]. In accordance with this algorithm, the thresholds and weights are iteratively adapted to minimize the error function E xð1Þ ;xð2Þ ;...;xðCÞ ;d ð1Þ ;d ð2Þ ;...;d ðCÞ ðw; hÞ ¼
C X M 1 X 2 ðyðiÞ d ðiÞ j Þ ; 2 i¼1 j¼1 j
(7) where w is the vector containing all weights in the neural net and h is the vector containing all thresholds of all ðiÞ neurons. The quantities yðiÞ ¼ ðyðiÞ 1 ; . . . ; yM Þ; i ¼ 1; . . . ; C; are the real output vectors of the neural net responding ðiÞ to the given set of the input vectors xðiÞ ¼ ðxðiÞ 1 ; . . . ; xN Þ ðiÞ (denoted as a learning set), d ðiÞ ¼ ðd 1ðiÞ ; . . . ; d M Þ are the required (ideal) output vectors corresponding to the input vectors xðiÞ and C represents the range (number) of the learning set. The chosen learning set contains 240 input vectors constructed using the P/D area ratios of 120 selected healthy OD images and of 120 selected glaucomatous OD images with the required neural net output binary values equal to the minimal, respectively,
ARTICLE IN PRESS 138
F. Pluha´cˇek, J. Pospı´ sˇ il / Optik 116 (2005) 133–143
maximal possible value 0, respectively, 1 of the neural net output. The considered iterative learning algorithm starts with the initial values w(0), hð0Þ: Simultaneously, the single components of the initial vectors w and h should be low [19,20], e.g. chosen at random from the interval (0.3, +0.3). In the lth step of the considered learning process, the values of weights and thresholds change from the actual values w(l) and hðlÞ to the new values w(l+1) and hðl þ 1Þ according to the relations wðl þ 1Þ ¼ wðlÞ Z1 gw þ a1 ½wðlÞ wðl 1Þ;
(8)
hðl þ 1Þ ¼ hðlÞ Z2 gy þ a2 ½hðlÞ hðl 1Þ;
(9)
where the symbols gw, respectively, gy represent the vectors of partial differentiations of the error function (7) with respect to the weights, respectively, thresholds, of the neurons in all layers of the neural net and Z1 ; Z2 ; a1 ; a2 are suitable learning parameters. We realized 50 learning iterations with empirically determined values Z1 ¼ Z2 ¼ 0:1; a1 ¼ a2 ¼ 0:005: The dependence of the error function (7) on the number of iterations of the utilized learning process is shown in Fig. 5. Simultaneously, the used values of w and h; also named the learned values, relate to the considered last iterative step. After the above-mentioned learning, the considered neural net was tested [19,20] using the test set of 240 input vectors constructed from the P/D area ratios of 120 selected healthy eye images, 120 selected glaucomatous eye images and the corresponding required output values that were not used in the learning set. According to the considered learning set mentioned above, the required neural net output was chosen equal to 0 for healthy eye images and 1 for glaucomatous eye images. During the test, the value of the error function (7) was computed for both the introduced test set and the used learning set. The obtained difference DE between the
Fig. 5. Dependence of the error function (7) on the number of iterations of the utilized learning process.
error function value for the test set and the learning set was acceptably small, DE 0:18: In view of the nonzero single standard deviations sT ; sI ; sII ; sIII and sIV ; the corresponding mean standard deviation sn 0:02 of the neural net outputs was numerically estimated on the basis of the common error-propagation rule [17].
5.3. Proposal of the neural net eye diagnostic classification and its analysis For the sought P/D area ratios of all investigated 240 healthy eye images and 240 glaucomatous eye images, the ð2Þ ð240Þ corresponding real neural net outputs yð1Þ h ; yh ; . . . ; y h ð1Þ ð2Þ ð240Þ and yg ; yg ; . . . ; yg were computed. The obtained histograms of the statistical distribution of neural net outputs yh in the case of healthy eye images and of the statistical distribution of neural net outputs yg in the case of the glaucomatous eye images are shown in Figs. 6 and 7. The eyes, investigated by the described neural net method, were sorted out into the three classes: lowglaucoma-risk class, glaucoma-suspected class and highglaucoma risk class. The low-glaucoma-risk class includes eyes with the corresponding neural net outputs below t1, the glaucoma-suspected class involves eyes with the corresponding neural net outputs within the interval ht1 ; t2 i and the high-glaucoma-risk class includes eyes with the corresponding neural net outputs above t2. At the same time, the parametric values t1 and t2, 0ot1 ot2 o1; are the suitable distinguishing limits. The chosen values t1 ¼ 0:13; t2 ¼ 0:63 are considered in this chapter. Thus the probability Ph/high, that the healthy eye is inconveniently classified to the high-glaucoma-risk class, is Ph=high ¼ 1 Dh ðt2 Þ ¼ 0:05
(10)
and the probability Pg/low, that the glaucomatous eye is inconveniently classified to the low-glaucoma-risk
Fig. 6. Obtained histogram of the statistical distributions of the neural net outputs yh in the case of healthy eye images.
ARTICLE IN PRESS F. Pluha´cˇek, J. Pospı´ sˇ il / Optik 116 (2005) 133–143
139
which expressions follow from the common errorpropagation rule of the explanation in [17]. The concrete values sh=high 0:01; sg=low 0:02; sh=suspected 0:04; sg=suspected 0:02 were numerically obtained by computation for the considered distinguishing values t1 ¼ 0:13; t2 ¼ 0:63; marked in Figs. 6 and 7.
6. Probability density curves method for diagnostic evaluation of the total P/DT area ratios 6.1. Probability densities estimation of the P/DT statistical distributions
Fig. 7. Obtained histogram of the statistical distributions of the neural net outputs yg in the case of glaucomatous eye images.
The probability densities ph and pg of the experimentally obtained statistical distributions of the x ¼ P=DT values of the above-mentioned 240 healthy and 240 glaucomatous OD images can be approximated by the experimentally ascertained relations
class, is
ph ðP=DT Þ ¼ N h exp½bh1 ðP=DT bh2 Þ3 ; P=DT 2 h0; 1i; (18)
Pg=low ¼ Dg ðt1 Þ ¼ 0:05:
(11)
The above-introduced functions Dh(t) and Dg(t) are experimentally obtained statistical distribution functions of the neural net outputs yh and yg. On the other hand, the probability Ph/suspected, respectively, Pg/suspected, for the classification of healthy, respectively, glaucomatous, eye to the glaucoma-suspected class can be computed by the relations Ph=suspected ¼ Dh ðt2 Þ Dh ðt1 Þ ¼ 0:45;
(12)
pg ðP=DT Þ ¼ N g exp½bg1 ðP=DT bg2 Þ2 ; P=DT 2 h0; 1i; (19) where Z
1
Nh ¼ 1
exp½bh1 ðx bh2 Þ3 dx;
(20)
exp½bg1 ðx bg2 Þ2 dx
(21)
0
Z
1
Ng ¼ 1 0
Pg=suspected ¼ Dg ðt2 Þ Dg ðt1 Þ ¼ 0:34:
(13)
With respect to the already mentioned non-zero mean standard deviation sn of the neural net outputs, the introduced probabilities Ph/high, Pg/low, Ph/suspected and Pg/suspected were determined with the relevant standard deviations @Ph=high (14) sn ; sh=high @t2 @Pg=low sg=low (15) sn ; @t1 " sh=suspected sn
@Ph=suspected @t1
2
#1=2 @Ph=suspected 2 þ ; @t2
are the real normalization constants and bh1 ; bh2 ; bg1 ; bg2 (shortly denoted by bh1;2 and bg1;2 ) represent the suitable real parameters. The functions (18) and (19) were found by histograms of P/DT values of healthy and glaucomatous eye images using the least-squares mathematical method described in [17,18]. For each histogram, the 10 class intervals were utilized. All computations were performed numerically with a sufficient accuracy. The obtained estimations of parameters bh1;2 ; bg1;2 and their standard deviations sh1 ; sh2 ; Table 1. Using least-squares method obtained estimations of parameters bh1;2 ; bg1;2 and their standard deviations sh1;2 ; sg1;2 Parameter
Estimated value
Standard deviation
Estimated value
bh1 bh2 bg1 bg2
7.0 102 11.4 102 5.2 101 24.3 102
sh1 sh2 sg1 sg2
0.6 102 0.4 102 0.4 101 0.4 102
(16) sg=suspected
"
#1=2 @Pg=suspected 2 @Pg=suspected 2 sn þ ; @t1 @t2 (17)
ARTICLE IN PRESS 140
F. Pluha´cˇek, J. Pospı´ sˇ il / Optik 116 (2005) 133–143
Fig. 8. Resultant histogram of the statistical distribution of the healthy eye P/DT images with the probability density curve ph(P/DT) of expression (18).
Fig. 9. Resultant histogram of the statistical distribution of the glaucomatous eye P/DT images with the probability density pg(P/DT) of expression (19).
sg1 ; sg2 (shortly expressed by sh1;2 and sg1;2 ) are introduced in Table 1. Both these probability densities (18) and (19) were also tested by the chi-squared test of fit [17,18] with the result that they did not reject for the significant levels below 22% in the case of ph and below 42% in the case of pg. The obtained resultant histograms of the healthy and glaucomatous eye images P/DT with the corresponding curves of the probability densities ph, pg are presented in Figs. 8 and 9.
6.2. Proposal of the probability density curves for the eye diagnostic classification and its analysis For the eye diagnostic classification, based on the described probability density curves method, the following three glaucoma risk classes of eyes were considered in the modified analogy with the previous neural net method: low-glaucoma-risk class, glaucomasuspected class and high-glaucoma-risk class. The
Fig. 10. Graphical dependences of the probability densities ph and pg on the P/DT area ratios under validity of the relations (18) and (19) for the in Table 1 introduced values of bh1;2 and bg1;2 :
low-glaucoma-risk class includes the eyes with the P/DT below t1, the glaucoma-suspected class involves the eyes with the P/DT from the interval ht1 ; t2 i and the high-glaucoma-risk class includes the eyes with the P/DT above t2. The distinguishing parametric limits t1 and t2, 0ot1 ot2 o1 are the suitably chosen values of P/DT, see Fig. 10. Under these conditions, the healthy eye is inconveniently classified to the high-glaucoma-risk class with the probability Z 1 Ph=high ¼ ph ðxÞ dx; (22) t2
while the glaucomatous eye is inconveniently classified to the low-glaucoma-risk class with the probability Z t1 Pg=low ¼ pg ðxÞ dx: (23) 0
With respect to the nonzero values of the standard deviations sT ; sh1;2 and sg1;2 ; the introduced probabilities Ph/high and Pg/low were determined with the relevant standard deviations sh=high ; sg=low of approximate definitions @Ph=high sT ; (24) sh=high @t2 @Pg=low sg=low sT ; (25) @t1 which follow from the common error-propagation rule. The performed numerical analysis showed that the influence of sh1;2 and sg1;2 on sh=high ; sg=low is negligible and therefore need not be considered. For the chosen values t1 ¼ 0:09; t2 ¼ 0:22; the equalities Ph=high ¼ 0:05; sh=high 0:02; Pg=low ¼ 0:05; sg=low 0:01 are valid.
ARTICLE IN PRESS F. Pluha´cˇek, J. Pospı´ sˇ il / Optik 116 (2005) 133–143
141
The remaining probability Ph/suspected for the classification of a healthy eye in the glaucoma-suspected class and the remaining probability Pg/suspected for the classification of a glaucomatous eye in the glaucomasuspected class now fulfil the relations Z t2 Ph=suspected ¼ ph ðxÞ dx; (26) t1
Z
t2
pg ðxÞ dx:
Pg=suspected ¼
(27)
t1
The corresponding relevant standard deviations sh=suspected ; sg=suspected for Ph/suspected and Pg/suspected can be computed using again the above-mentioned common error-propagation rule, which gives the equations "
#1=2 @Ph=suspected 2 @Ph=suspected 2 sh=suspected sT þ ; @t1 @t2 (28)
sg=suspected
"
#1=2 @Pg=suspected 2 @Pg=suspected 2 sT þ : @t1 @t2 (29)
In our case for the considered equalities t1 ¼ 0:09 and t2 ¼ 0:22; the values Ph=suspected ¼ 0:45; sh=suspected 0:04; Pg=suspected ¼ 0:34 and sg=suspected 0:03 were obtained.
6.3. Glaucoma prevalence probability in view of the total P/DT area ratios and the subject ages For the total probability P of the occurrence of glaucomatous eyes in a population group of eyes with the given P/DT value and its standard deviation sT (in short called the glaucoma prevalence probability), the corresponding formula P¼
P p pg Pp pg þ ð1 Pp Þph
(30)
was derived on the basis of the conditional probability theory [18]. In our case, the used probability densities ph and pg are given by the relations (18) and (19). The introduced parameter Pp with the standard deviation sPp means the probability of glaucoma occurrence in the whole population group (without reference to the P/DT values) of the given age range. In general, Pp is markedly influenced by the race of the investigated subjects and differs for different age ranges, see [1,21]. The relevant resultant standard deviation sP for the total probability P was estimated according to the common errorpropagation rule by the suitable relation
Fig. 11. Graphs of the quantities (30) and (31) versus the P/DT area ratios constructed for the examined subject age ranges below 65 years and above 75 years in the case of the parametric values of bh1;2 ; bg1;2 and sh1;2 ; sg1;2 introduced in Table 1.
" sP
@P @ðP=DT Þ
2
s2T
@P þ @Pp
2
#1=2 s2Pp
;
(31)
where the influence of sh1;2 and sg1;2 was acceptably neglected. Relations (30) and (31) are plotted as functions of P/DT in Fig. 11 for two different age ranges of the white population (below 65 years and above 75 years) and for in Table 1 introduced estimated values of bh1;2 ; bg1;2 ; sh1;2 and sg1;2 : In the case of age below 65 years, the values Pp ¼ 0:01; sPp ¼ 0:005 were obtained and for age above 75 years, the equalities Pp ¼ 0:03; sPp ¼ 0:01 are valid, see [21]. It is evident from the graphs in Fig. 11 that the glaucoma prevalence probability P increases with the increased P/DT values in both considered cases. Moreover, the resultant standard deviations sP show the maximum and their values are inconveniently large, mainly for the interval of the maximal slope of P.
7. Comparison of the presented methods and the obtained results The described quantile curves method, based on the utility of the quantile standpoint for diagnostic evaluation of the healthy and glaucomatous human eyes, enables advantageously the convenient direct objective determination of the corresponding quantile curves and then also of the number of the representative OD sectors with undesirably increased values of the P/D area ratios. The disadvantage of this method is that the final glaucoma diagnosis must be often accompanied with a subjective decision without any possibility to compute directly the resultant diagnostic probability. On the other hand, the also presented neural net and probability density curves methods enable
ARTICLE IN PRESS F. Pluha´cˇek, J. Pospı´ sˇ il / Optik 116 (2005) 133–143
142
advantageously an objective final automatic classification of the investigated eyes in three exactly defined glaucoma risk classes. Because probability values Ph/high or Pg/low for classification of healthy or glaucomatous eye in the high-glaucoma or low-glaucoma risk class are conveniently small enough for both later-mentioned glaucoma diagnostic methods (Ph=high ¼ Pg=low ¼ 0:05), the eyes from the high-glaucoma risk class, respectively, low-glaucoma risk class, can be reliably diagnosed as really glaucomatous, or healthy. In such cases, the probability of incorrect diagnosis is namely also equal to Ph/high, or Pg/low, i.e. it is small. The remaining probabilities Ph/suspected and Pg/suspected for classification of healthy or glaucomatous eye in the glaucomasuspected class are both larger (Ph=suspected ¼ 0:45 and Pg=suspected ¼ 0:34). Thus a reliable diagnosis of the eye must be carried out differently, by means of a characteristic different quantity. According to Machala and Pospı´ sˇ il [22], the suitable reliability factor of possible definition r¼
j¯th ¯tg j ½ðs2h
þ s2g Þ=21=2
(32)
can be applied. The corresponding characteristic parameters ¯th ; ¯tg and sh ; sg for healthy and glaucomatous eyes are the mean values and standard deviations relating to the obtained statistical distribution functions Dh(t) and Dg(t) for the outputs t ¼ yh and yg in the case of the neural net method, or the mean values and standard deviations of the obtained distributions of the probability densities ph(t) and pg(t) for the values t ¼ P=DT in the case of the probability density curves method. From the definition (32) it is evident that the value of r; representing in fact the normalized distance between ¯th and ¯tg ; increases with increased difference j¯th ¯tg j and decreases with increased standard deviations sh ; sg : A higher value of r corresponds to a higher diagnostic reliability of the used method and vice versa. Values r ¼ 1:82 and 1.72 were obtained for the neural net method and the probability density curves method. Because these values differ only little for both considered methods, they have an approximately identical reliability. In addition, the probability density curves method enables the computation of the glaucoma prevalence probability for the given total P/DT area ratios and the subject age ranges. With respect to these facts and to the complicated structure of the utilized neural net, the probability density curves method is mostly more useful in practical applications.
8. Conclusion The contribution of the present article consists in the proposal, verification and comparison of three new
computer image analysis statistical methods for experimental diagnostic evaluation of colour glaucomatous changes within the optic disc (OD) of a human eye retina. The search for such changes by means of a suitable computer-connected digital fundus camera is based on the preliminary determination of the total representative OD areas and their relevant relative pallor area sizes (P/D), which are represented by 36dimensional characteristic vectors of components xn ¼ P=D1 ; P=D2 ; . . . ; P=D36 or, respectively, by five-dimensional characteristic vectors of components xk ¼ P=DT ; P/DI, P/DII, P/DIII, P/DVI. The suitable established components relate to the elementary P/D area ratios for 36 equally spaced OD sectors or, respectively, to the combination of the total OD area (P/DT) and the partial quadrants numbered I, II, III and IV in the total OD area for the sufficiently large representative sets of the normal healthy and the glaucomatous eyes. As the first is presented the so-called the quatnile curves method, based on the knowledge of the vector components P/Dn, n ¼ 1; 2; . . . ; 36; and their experimentally obtained distribution functions Dhn(P/Dn) and Dgn(P/Dn) for the sets of 240 normal healthy and 240 glaucomatous eyes and for the special parametric values ðxÞ qðxÞ hn and qgn of P/Dn relating to the suitable chosen normal and glaucomatous quantile x: The resultant main dependences of these well-arranged values in the descending order on their sequence in this ordering represent the x% normal and glaucomatous quantile curves for the successive and effective evaluation purposes also in comparison with the suitably introduced one-eye normal and glaucomatous curves. The second investigated method, specified as the neural net method, utilizes the modified artificial computer-simulated neural net point of view with the corresponding back-propagation learning algorithm for diagnostic evaluation of the 240 normal healthy and 240 glaucomatous eyes on the basis of their fivedimensional characteristic input vectors of type x ¼ ðP=DT ; P=DI ; P=DII ; P=DIII ; P=DIV Þ: For the searched P/D area ratios of all investigated healthy and glaucomatous eye images, the corresponding resultant neural ð2Þ ð240Þ ð2Þ ð240Þ net outputs yð1Þ and yð1Þ g ; yg ; . . . ; yg h ; yh ; . . . ; yh were computed and their histograms were constructed. On the basis of the resultant histograms, the investigated eyes were classified into three glaucoma risk classes and diagnosed by means of the significant error probabilities and standard deviations. The last considered method, i.e. the probability density curves method, relates to the obtained statistical distributions of the P/DT values for the also considered 240 healthy and 240 glaucomatous OD images. It is based on the exploitation of the experimentally ascertained relations for the probability density curves ph(P/DT) and pg(P/DT) relating to the mentioned sets of the investigated normal healthy and glaucomatous
ARTICLE IN PRESS F. Pluha´cˇek, J. Pospı´ sˇ il / Optik 116 (2005) 133–143
eyes and to a suitably chosen parameters. These curves were found by the corresponding histograms of the P/ DT values using the least-squares mathematical method. On the basis of these curves, the investigated eyes were classified into three glaucoma risk classes, on analogy with the previous neural net method, but by means of the significant error probabilities and standard deviations of different concrete definitions. Moreover, also the total glaucoma prevalence probability and the corresponding resultant standard deviation were estimated for different age ranges of the investigated subjects. Both the last two automatic and objective statistical diagnostic methods can replace effectively a subjective evaluation of glaucomatous changes within the representative OD. The performed analysis shows that the reliabilities of these methods are approximately identical. In comparison with other classical non-statistical glaucoma diagnostic methods, based directly on the numerical representation of the OD glaucomatous changes (e.g. the stereoscopic images computer evaluation method or the laser scanning tomography method), the designed statistical methods require only the computer-connected digital fundus camera arrangement without other basic devices. They enable also a larger diversity of the current ophthalmologic equipment without any new expensive investment.
Acknowledgement All the investigated optic disc images were kindly supplied by the physician MUDr. Toma´sˇ Kubeˇna. The authors would like to express their gratitude to him.
References [1] B. Shields, Textbook of Glaucoma, Williams & Wilkins, Baltimore, 1992. [2] H. Kraus, J. Boguszakova´, P. Diblı´ k, D. Dotrˇ elova´, M. Filipec, J. Hycl, I. Karel, I. Kocour, H. Kraus, J. Krˇ epelkova´, P. Kuchynka, J. Otradovec, J. Palecˇkova´, J. Peregrin, J. R˘ehu˚rˇ ek, J. Sveˇra´k, A. Sˇedivy`, Compendium of the Eye Medicine, Grada, Praha, 1997 (in Czech). [3] A.P. Nesterov, Primary Glaucoma, Avicenum, Praha, 1991 (in Czech). [4] S. R˘eha´k, R. Knobloch, O. Riebel, F. Vrabec, Medicine of the Eye, Avicenum, Praha, 1980 (in Czech). [5] A. Nover, The Ocular Fundus. Methods of Examination and Typical Findings, Schattauer, Stuttgart, 1987.
143
[6] H. Sauter, W. Straub, H. Rossmann, Atlas of the Ocular Fundus. Photographs of Typical Changes in Ocular and Systemic Disease, Urban & Schwarzenberg, Mu¨nchen, 1997. [7] T.N. Cornsweet, S. Hersh, J.C. Humphries, R.J. Beesmer, D.W. Cornsweet, Quantification of the shape and color of the optic nerve head, In: G.M. Breinin, I.M. Siegel (Eds.), Advances in Diagnostic Visual Optics, Springer, New York, 1983, pp. 141–149. [8] R. Varma, G.L. Spaeth, The par is 2000: a new system for retinal digital image analysis, Ophthal. Surg. 19 (1988) 183–192. [9] T.N. Cornsweet, S. Hersh, Ocular-fundus analyzer. US Patent No. 4715703, US Government, Washington, 1987. [10] F. Pluha´cˇek, J. Pospı´ sˇ il, Present methods of the computer-image analysis of a human retina with regard to diagnostics of glaucoma, Fine Mech. Opt. 46 (2001) 326–334 (in Czech). [11] T. Shioiri, Ophthalmologic image processor, US Patent No. 6104828, US Government, Washington, 2000. [12] F.S. Mikelberg, K. Wijsman, M. Schulzer, Reproducibility of topographic parameters obtained with the heidelberg retina tomograph, J. Glaucoma 2 (1993) 101–103. [13] G. Zinser, U. Harbarth, H. Schro¨der, Formation and analysis of three-dimensional data with the laser tomographic scanner (LST), In: J.E. Nasemann, R.O.W. Burk (Eds.), Scanning Laser Ophthalmoscopy and Tomography, Quintessenz Verlag, Mu¨nchen, 1990, pp. 243–252. [14] K.U. Bartz-Schmidt, A. Sengersdorf, P. Esser, P. Walter, R.D. Hilgers, G.K. Krieglstein, The cumulative normalised rim/disc area ratio curve, Greafe’s Arch. Clin. Exp. Ophthalmol. 234 (1995) 227–231. [15] J. Caprioli, Reproducibility of optic disc measurements with computized analysis of stereoscopic video images, Arch. Ophthalmol. 104 (1986) 1035–1039. [16] T. Kubeˇna, P. Cˇernosˇ ek, J. Mayer, Technique of digital planimetry of optic disc, Czech Slovak Ophthalmol. 56 (2000) 170–175 (in Czech). [17] L. Kuba´cˇek, L. Kuba´cˇkova´, Statistics and Metrology, Palacky` University Press, Olomouc, 2000 (in Czech). [18] J. Sˇkra´sˇ ek, Z. Tichy`, Principles of Applied Mathematics III, SNTL, Praha, 1990 (in Czech). [19] B. Mu¨ller, J. Reinhardt, M.T. Strickland, Neural Networks: An Introduction, Springer, Berlin, 1995. [20] J.G. Taylor, Neural Networks and their Applications, Wiley, Chichester, 1996. [21] B.E. Klein, R. Klein, W.E. Sponsel, T. Franke, L.B. Cantor, J. Martone, M.J. Menage, Prevalence of glaucoma: the beaver dam eye study, Ophthalmol. 99 (1992) 1499–1504. [22] L. Machala, J. Pospı´ sˇ il, Proposal and verification of two methods for evaluation of the human iris video-camera images, Optik 112 (2001) 335–340.