Automatic method of analysis of OCT images in assessing the severity degree of glaucoma and the visual field loss

Introduction In many practical aspects of ophthalmology, it is necessary to assess the severity degree of glaucoma in cases where, for various reasons, it is impossible to perform a visual field test - static perimetry. These are cases in which the visual field test result is not reliable, e.g. advanced AMD (Age-related Macular Degeneration). In these cases, there is a need to determine the severity of glaucoma, mainly on the basis of optic nerve head (ONH) and retinal nerve fibre layer (RNFL) structure. OCT is one of the diagnostic methods capable of analysing changes in both, ONH and RNFL in glaucoma. Material and method OCT images of the eye fundus of 55 patients (110 eyes) were obtained from the SOCT Copernicus (Optopol Tech. SA, Zawiercie, Poland). The authors proposed a new method for automatic determination of the RNFL (retinal nerve fibre layer) and other parameters using: mathematical morphology and profiled segmentation based on morphometric information of the eye fundus. A quantitative ratio of the quality of the optic disk and RNFL – BGA (biomorphological glaucoma advancement) was also proposed. The obtained results were compared with the results obtained from a static perimeter. Results Correlations between the known parameters of the optic disk as well as those suggested by the authors and the results obtained from static perimetry were calculated. The result of correlation with the static perimetry was 0.78 for the existing methods of image analysis and 0.86 for the proposed method. Practical usefulness of the proposed ratio BGA and the impact of the three most important features on the result were assessed. The following results of correlation for the three proposed classes were obtained: cup/disk diameter 0.84, disk diameter 0.97 and the RNFL 1.0. Thus, analysis of the supposed visual field result in the case of glaucoma is possible based only on OCT images of the eye fundus. Conclusions The calculations and analyses performed with the proposed algorithm and BGA ratio confirm that it is possible to calculate supposed mean defect (MD) of the visual field test based on OCT images of the eye fundus.


Introduction
Analysis of OCT images of the eye fundus is a technique known for many years which enables to acquire large amounts of information useful for medical diagnosis. Today, there are many known profiled algorithms for the analysis of both successive layers of the eye fundus [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20] as well as the choroid layer [21]. Depending on the approach and practical implementation, these algorithms operate in the area of texture analysis [1,2], algorithms based on segmentation with Canny method [4], SVM [5], mathematical morphology [6], active contour [7,8], fuzzy algorithms [9], Markov model [10], robust segmentation [11], random contour analysis [12] and others [13][14][15][16][17][18][19][20]. Implementation of these algorithms most commonly occurs in the C, C# programming environment if they are used in commercially available applications, or, for example, in Matlab if they are at the stage of testing and modification [6]. The analysis time of a single OCT image of the fundus is several tens of milliseconds to a few seconds with the most computationally intensive applications [1,17]. In practice, the biggest difficulty is large interindividual variability and a possible high degree of retinal pathology. This translates into problems with the designation of successive layers of the eye fundus (RNFL-retinal nerve fibre layer, RPE-retinal pigment epithelium, and others) in the following cases ( Figure 1): The layer (u 1 , u 3 ) is not continuous - Figure 1a  The following issues are shown in sequence: a) the problem with the gap between u 1 and u 3 and the impact of the additional new layer u 2 , b) the problem with combining layers u 2 and u 3 with the layer u 1 for 3 pixels (1,4), (1,5) and (1,6), c) the difficulty in deciding which layer u 1 or u 2 should be combined with the layer u 3 , d) the difficulty in deciding whether the layer u 2 is the missing piece connecting layers u 1 and u 3 .
The mentioned problems can be handled with in a limited way using artificial intelligence (knowledge gained from other images) or anthropometric evidence. In any case, image analysis and the proposed approach to solving such cases are difficult.
From a medical viewpoint, the RNFL thickness is essential in the case of the optic disk analysis. On the basis of its average thickness and further measurements, other parameters such as a cup/disk area ratio or a cup/disk vertical diameter are obtained. These parameters have an influence on the diagnosis of an ophthalmologist as well as the applied methods of treatment.
In the case of glaucoma, it is essential to test the visual field. Its measurement is usually carried out using a static perimeter. A reliable parameter obtained from perimetric data is the average depth of a defect. In many cases, for example, in a number of different additional diseases, it is impossible to obtain this type of measurement or the results obtained are not reliable. Then, it is necessary to search for other methods which will enable to determine the supposed visual field results. One such method is the analysis of OCT images of the eye fundus proposed by the authors. It allows not only to determine the severity degree of glaucoma but also to determine the supposed MD value of visual field test. It is further shown which selected parameters in OCT images correlate, to the greatest extent, with the results obtained from perimetry. Moreover, the new algorithm proposed by the authors and the quantitative ratio of the optic disk condition are described.

Material
OCT images of the eye fundus of 55 glaucoma patients (110 eyes) were obtained from the SOCT Copernicus (Optopol Tech. SA, Zawiercie, Poland). Patients ranged in age from 37 to 88 years. The obtained images (B-scans) had a resolution of M × N = 1000 × 600 pixels (where Mnumber of rows, Nnumber of columns) at a colour resolution of 8 bits/pixel. For each patient I = 99 B-scans were obtained allowing for full 3D reconstruction of the eye fundus image. At the same time, 110 results were acquired from the static perimeter OCTOPUS 300Series V 6.07c. All data were analysed in a source format. The presented analysis was carried out retrospectively on the basis of data obtained during typical medical examinations. The data were anonymised and the patients were examined in accordance with the Declaration of Helsinki.

Methods
The method of analysis and processing of OCT images proposed by the authors was preceded by the analysis of the results obtained from conventional methods. This analysis indicates features which should be taken into account when constructing the target image analysis algorithm.

Selection of the features of the image
As stated in the previous section, a visual field test was performed on 55 glaucoma patients. The features of the eye fundus were measured using a Copernicus tomograph. The results obtained in the software available for these devices are shown in Figure 2 and Figure 3. The features obtained on this basis are presented in Table 1. Examples of values for each feature for the first 10 patients are shown in Table 2. Each of these characteristics will be further, for simplicity, denoted by the symbol "w". In this case, the features from w(1) to w(28) were acquired from the tomograph, and w(29) and w(30) from the perimeter. The parameters from w(1) to w(30) were further subjected to exploratory factor analysis in order to detect their cross correlation and to remove those features which are statistically insignificant. For this purpose, the following items were designated in succession: correlation of features w, eigenvalues of the correlation matrix, the number of classes using the Kaiser's and Cattell's criteria, rotation of the coordinate system, correlation of individual features with the newly created classes. The results (for the whole group of patients) are shown in Table 3 whereas the correlation chart of various features with three classes are shown in Figure 4. The analysis of Table 3 and Figure 4 suggests several conclusions useful for further analysis: Class 1 significantly (above 0.5) correlates with the features w(16), w(13), w(6), w(5), w(8), w(3), w(15) and w (17) and is in a negative correlation with the features w(19), w(7) and w (4). It means that the cup/disk diameter and the cup diameter are a representative of the class 1. The ratio of the cup vertical/horizontal and the area and volume of the rim are in the inverse correlation.
Class 3 correlates with the features w(10), w(2), w(12), w(11) and w(14) which are responsible for the parameters of the disk diameter in the vertical and horizontal axes.
Analysis of this division into classes, their correlation with the results obtained from perimetry and interpretation of individual features allow to specify the features of OCT images necessary for further analysis. The particularly important features of OCT image are: disk diameter and its relation to cup as well as the RNFL parameters in the area. This information is the basis for the design of a profiled image analysis algorithm, correlation analysis with the results of perimetry and for suggesting a new quantitative ratio of the optic disk condition.

Figure 2
The result for one patient obtained from the Copernicus tomograph. The following images can be seen in the corners: fundus reconstruction, a deviation from the normal map, the RNFL thickness map and the graph RNFL TSNIT. In the central part at the bottom there is a sample B-scan of the optic disk. The central part is dedicated to the ratios calculated by the software.
Pre-processing Figure 3 The result for one patient obtained from the static perimeter. From the top: data on patient and study, a colour picture of the quality of vision in the 30 O , numerical data of quality of vision in various fields expressed in decibels: values, comparisons, corrected comparisons, probability and corrected probability. In the lower right corner there are two most important parameters, namely MD (Mean Defect) and LV (Loss Variance). The analysis did not take into account LV value, because there was no significant correlation between LV and any of the OCT parameters.
of A-scans and B-scans, date of test, the patient data and others. The sequence of images L GRAY (m,n,i) obtained on this basis is further subjected to filtration with a median filter whose mask h 1 is sized M h1 × N h1 × I h1 = 3 × 3 × 3 pixels. The mask size was chosen arbitrarily taking into account the size of distortions and artefacts entering the optical path. The filtered image L M (m,n,i) ( Figure 5) is subjected to further preliminary transformations. These include determination of the RNFL.
Processing -RPE, ILM, RNFL, cup, disk In order to determine the RNFL, the RPE layer (retinal pigment epithelium) is pre-determined automatically. The RPE layer is the most characteristic and the brightest layer in an OCT image, therefore, it is used to determine the target RNFL. Additionally, the RPE is used to determine the disk boundary. For this purpose, every n-th column of the image L M , each A-scan, is analysed. Analogously to the method described in [21,22], the position of the maximum brightness for each column is predetermined, i.e.: where: m,ncoordinates of the rows and columns of the matrix m∈(1,M) and n∈(1,N).
The formula (1) can be directly applied only when for all the analysed rows of each column there is only one maximum brightness value. As shown in [6], in practice it occurs in about 80% of the cases analysed at a resolution of 8 bits per pixel. On this where: meanmean value, p mthreshold equal 2/3. The value of 2/3 (p m ) was chosen arbitrarily and is associated with the maximum acceptable shift of the RPE layer in the area of the optic disk. The other element, namely the ILM layer area, was calculated as the area from the first pixel of each row to L RPE (n,i). It is very easy to perform binarization in this area using Otsu method [23], thus setting the binarization threshold p r . The resulting binary image L B (m,n,i) is: Based on the binary image L B (m,n,i), the edge of the ILM layer was defined as: L ILM n; i ð Þ ¼ min The waveform of L ILM (n,i) and the RPE boundary, i.e. L RPE (n,i), are the basis for further analysis of the RNFL layer and cup area. The area between the layers ILM and RPE is further subjected to an operation which involves determination of the edge The protection (m > L ILM (n,i))∧(m < L RPE (n,i)) is additional and prevents detection of the RNFL on the edge of the analysed object with the background. The resulting edge has a number of disadvantages. The biggest one is the lack of continuity caused by both occurring shadows as well as insufficient contrast in the place of its detection. For such cases, the method of active edge correction was suggested. This method involves adjustment of the position of the layer L RPE (n,i). For each point of the layer L RPE (n,i), the mean brightness value above and below L RPE (n,i) in the areas L u and L d is determined - Figure 6. The new corrected position of the layer L RPE * (n,i) is determined as the maximum difference between the average values in the areas L u (Δud,n,i) and L d (Δud,n,i). These areas are shifted in the column axis in the range Δud = ±med(|L ILM (n,i)-L RPE (n,i)|). The adjusted position of the RNFL is therefore calculated as: The resulting new position of the layer L RNFL * (n,i) is the basis for further analysis. On the basis of the positions of the ILM and RPE layers, i.e.: L ILM (n,i), L RPE (n,i), the cup area L cup (n,i) is determined. This area is calculated as the intersection of the boundary of L disk (n,i) at the offset height equal to 150 μm with the layer L ILM (n,i).
As evident from the exploratory factor analysis presented in the previous section, apart from the RNFL thickness, also limitation of the analysis area -ROI (region of interest) is vital. The ROI ranges from 2•r disk to 3•r disk in an OCT image and in order to determine it, it is necessary to determine the position of coordinates of the centre (n d0 ,i d0 ), i.e.: The value of r disk is the average radius calculated from the area of a circle: The calculation of the centre of gravity is followed by the conversion of the ROI to the polar coordinate system, i.e.: The function atan2 returns θ the same size as n and i containing the element-byelement, four-quadrant inverse tangent (arctangent) of the real parts of n and i. The obtained results, the ROI in the polar coordinate system L RNFL * (θ,ρ), are shown in Figure 7. The proposed algorithm is shown globally in a block diagram form in Figure 8. On the basis of the presented ROI, the following results (features) were obtained: the disk diameter and its relation to the cup and RNFL parameters. Their thorough analysis and correlation with the results from perimetry are presented in the next section.

Results
The obtained result L RNFL * (θ,ρ), the ROI in the polar coordinate system (Figure 7), is the basis for further analysis. Based on the presented exploratory factor analysis and the proposed method of image analysis and processing, the features representing different classes ( Figure 4, Table 3), associated with L RNFL * (θ,ρ), were formulated in the following way: w g (1)the relative difference in the positions of centres of gravity calculated as: where: n d0 , i d0coordinates of the centre of gravity of the disk, n c0 , i c0coordinates of the centre of gravity of the cup, w g (2)relative, inverse, minimum distance between the boundaries of the cup and disk i.e.: where: r minminimum distance between the boundaries of the cup and disk. w g (3) -RNFL thickness deviation in the ROI (L RNFL ROI (θ)) for particular areas: temporal (TM), superior (SU), nasal (NA), interior (IN) from the standard thickness changes: where: L RNFL WZ (θ)the reference value of the RNFL thickness (the graph highlighted in green in Figure 2).
The mentioned features are components of the proposed biomorphological glaucoma advancement ratio (BGA) which is defined as: Figure 8 Block diagram of the algorithm proposed by the authors. The patient undergoes a fundus examination with the use of a tomograph and their visual field is measured with a perimeter. Then, the significance and correlation of the features obtained are calculated, thereby forming three classes. They are the basis for proposing a new analysis algorithm. After automatic determination of the layers RPE, ILM, RNFL and the cup and disk, the features w g (1), w g (2) and w g (3) are determined. In the final stage, the proposed BGA ratio and correlation with results obtained from perimetry are calculated.
After value substitution and simplification: Figure 9 presents a list of various types of cup arrangements relative to the disk and the impact of their sizes on the values of the features w g (1), w g (2) and w g (3) as well as w BGA . In addition, for comparison, the values of the DDLS (Disk Damage Likelihood Scale) were also shown [24]. A change in the position of the centre of gravity of the cup relative to the disk is expressed by the feature w g (1). The closer the centres of gravity of the cup are to each other, the closer the feature w g (1) is to zero. If the centres of gravity of the cup and disk are far apart, the extreme position of the cup is on the edge of the disk, then the value of w g (1) is close to one. In turn, the feature w g (2) is a measure of centricity of the distribution of cup edges relative to the disk. When the edge of Figure 9 The impact of cup arrangement with respect to the disk on the values w g and w BGA . A summary of various typical cup arrangements relative to the disk and the effect of their size on the features w g (1), w g (2) and w g (3) as well as w BGA and, in addition, the values of the DDLS, have been shown. A change in the position of the cup centre of gravity relative to the disk is expressed by the feature w g (1). In turn, the feature w g (2) is a measure of the centricity of the cup edge distribution with respect to the disk. w g (3) is the mean value of the RNFL thickness deviation from the standard. the cup is substantially spaced from the edge of the disk, the value of w g (2) is close to zero. When the edges of the cup and the disk overlap in any area, then the value of w g (2) is close to one. The situation is similar in the case of advanced glaucoma where r min is close to zero and w g (2) is close to one. A change in the value of w g (3) is related to a deviation with respect to the standard waveform of the RNFL thickness ( Figure 2). The greater this difference is, the larger the value of w g (3) is, and vice versa. The features w g (1), w g (2) and w g (3) are summed up and thus determine the ratio w BGA. The ratio w BGA reaches values close to one if all three features w g indicate pathology. In the case of the right features w g , the ratio w BGA is equal to zero.
Based on the proposed ratio w BGA , its correlation with the results of perimetry was determined. In addition, it was compared with the results obtained only on the basis of known software (features w(1) to w(29)). The following results were obtained: Correlation of classes 1, 2, 3 (obtained from features w(1) to w(29)) with the results obtained from perimetry is 0.78, Correlation of the proposed ratio w BGA referred to the described method of image analysis and processing with the results of perimetry is 0.86 - Figure 10. The function describing the dependence of w(29) (denoting the MD value from perimetry) on w BGA takes the following form: The difference in the correlation results stems from the profiled method of image analysis and processing and the measurement methodology of the described issue. Figure 10 Graph of the dependence of w(29) as a function of changes in the ratio w BGA . The graph shows in blue the results obtained for individual patients. The result of linear regression is marked in green. The function describing the dependence of w(29) on w BGA takes the following form: w(29) = 24⋅w BGA -9.34. On this basis, when the ratio w BGA is calculated, it is possible to calculate the expected result of perimetry.